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1 . 0  INTRODUCTION 
a.  General 

Progress  in  the  development  of  practical  structural  optimization 
technology  has  been  consistently  characterized  by  a  series  of  major  advances, 
followed  by  periods  of  consolidation  and  even  retrenchment.  Design  or 
optimization  is  fundamentally  a  more  sophisticated  concept  than  structural 
analysis  and  it  is  therefore  only  to  be  expected  that  the  derivation  and  develop¬ 
ment  of  new  techniques  will  be  a  much  slower  process  than,  say,  the  develop¬ 
ment  of  finite  element  analysis  methods. 

After  the  major  concentration  of  attention  in  the  60' s  on  mathematical 
programming  techniques  reached  what  was  effectively  a  computational 
stalemate,  the  development  of  optimality  criteria  methods^  in  the  early  70's 
appeared  to  offer  major  improvements  over  any  other  optimization  methods 
currently  in  vogue.  Certain  deficiencies  in  optimality  criteria  methods  were 
recognized  but  the  computational  advantages  of  some  approximate  methods 
were  sufficient  to  warrant  the  development  of  large  scale  computer  programs, 
such  as  OPTIM(2,  3)  and  ASOP(4). 

Following  the  initial  rush  of  development  has  been  a  larger  period  of 
study,  during  which  the  deeper  implications  and  above-mentioned  deficiencies 
in  the  approximate  optimality  criteria  methods  have  been  evaluated.  Although 
a  number  of  differing  approaches  to  redesign  strategies  have  been  proposed 
by  a  number  of  authors,  it  has  been  demonstrated  ^  5  ^  that  the  vast  majority 
of  these  approaches  are  really  variations  of  a  common  theme.  The  essential 


differences  being  merely  resident  in  the  selection  of  arbitrary  coefficients 
and  other  parameters  which  might  appear  to  influence  convergence 
characteristics. 

In  spite  of  this  work,  it  appears  that  little  real  headway  has 
been  made  in  solving  some  of  the  fundamental  problems  in  optimization. 

The  use  of  a  Lagrangian  formulation  is  well  established  as  a  starting  point 
for  optimization  but  there  are  still  many  outstanding  obstacles  between  the 
initial  concept  and  workable  optimization  tools.  For  single  displacement 
constraints,  there  is  no  problem.  For  multiple  displacement  constraints,  the 
governing  optimality  criteria  can  be  established,  but  the  solution  of  these 
nonlinear  equations  with  unknown  constraint  population  is  by  no  means  a 
solved  problem.  For  stress  constraints,  even  the  establishment  of  valid 
optimality  criteria  has  not  hitherto  been  successfully  performed.  In  summary, 
it  may  be  stated  that  there  are  no  adequate  operational  methods  for  optimizing 
structures  subject  to  multiple  constraints  on  strength  (stress),  stiffness 
(displacement)  and  fabricational  limits  (minimum  member  sizes).  Other  types 
of  constraints -stability,  dynamic  response,  flutter  and  aeroelasticity  are 
of  prime  importance,  but  they  are  generally  and  generically  related  to  the 
three  basic  constraints.  Hence,  the  approach  is  to  solve  the  problems 
associated  with  the  basic  constraints  and  the  more  sophisticated  constraints 
will  follow. 

For  stress  constraints,  the  classic,  time-honored  and  erroneous 
fully  stressed  design  (FSD)  method  based  upon  the  stress-ratio  technique 
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ia  most  widely  used.  FSD,  while  inaccurate,  has  the  overwhelming  merit 
of  being  extremely  simple  and  economical  to  use.  For  the  treatment  of 
displacement  constraints  no  equivalent  existed  in  the  60' s  and  only  rigorous 
(and  costly)  mathematical  programming  methods  could  be  used. 

The  development  of  the  methods  now  generally  labeled  as  optimality 
criteria  was  the  result  of  attempts  to  fill  this  need  for  an  approach  to 
displacement  constraints.  The  optimality  criteria  methods  are  iterative 
(as  are  mathematical  programming  methods)  but  the  number  of  iterations 
required  for  convergence  appeared  to  be  largely  independent  of  the  number 
of  variables  (the  downfall  of  pure  math  programming). 

A  considerable  degree  of  success  was  achieved,  particularly  using 
the  envelope  method,  but  the  limitations  of  the  approach  were  recognized. 

The  solution  is  exact  for  a  single  equality  (displacement)  constraint  only. 

For  multiple  constraints, the  procedure  for  active/passive  partitioning  of 
members  is  approximate  and  there  is  no  capability  for  treating  problems 
involving  mixed  equality  and  inequality  constraints.  In  addition,  the 
incorporation  of  stress  and  fabricational  constraints  is  effectively  based 
upon  the  FSD  method. 

Work  has  been  carried  on  by  a  number  of  research  teams  in  attempts 
to  (a)  solve  the  problem  of  multiple  displacement  constraints  in  a  rigorous 
(or  even  non-rigorous  but  practical)  manner  -  but  with  little  real  success  and 
(b)  improve  upon  stress-ratio  FSD  for  strength  constraints.  For  the  latter, 
concepts  such  as  strain  energy  density,  adaptive  steps,  etc.  ,  have  been 
investigated  but  none  of  the  approaches  have  been  totally  effective. 


The  finite  element  displacement  method  has  dominated  the  field  of 

analysis  for  many  years  and  has  been  used  almost  exclusively  in  optimization 

work.  The  force  method  has  largely  been  ignored  in  analysis  and  consequently 

has  not  been  considered  of  general  or  particular  relevance  to  optimization. 

Some  recent  use  has  been  made  of  the  force  method^'  but  this  was 

primarily  to  reduce  the  computational  effort  associated  with  iterative  analysis. 

(o) 

At  Bell,  the  role  of  the  force  method  was  examined'  to  determine  the 
possibility  of  a  fundamental  integration  of  the  analysis  and  redesign  philosophy. 
The  preliminary  results  did  indicate  a  potential  for  the  rigorous  incorporation 
of  stress  constraints.  While  extremely  simple,  the  approach  of  using  the 
force  method  concept  to  overcome  the  difficulties  associated  with  strength 
optimization  of  redundant  systems  appears  to  be  entirely  novel.  If  the 
internal  redundant  forces  are  considered  to  be  self-equilibriating  external 
forces,  the  structure  considered  becomes  effectively  statically  determinate. 
The  optimality  criteria  for  such  a  determinate  structure  can  then  be  derived 
rigorously  for  either  displacement  or  stress  constraints.  The  compatibility 
conditions,  associated  with  the  redundancies,  now  become  zero-valued 
quasi-external  displacement  constraints.  Additional  optimality  criteria  are 
derived  to  permit  the  consideration  of  the  redundant  forces  as  subsidiary 
variables. 

The  problem  of  stress  constraints  is  thus  transformed  into  one 
involving  displacement  constraints.  This  general  concept  has  also  been 
proposed  for  bar-type  structures^  ^  \  whereby  the  maximum  allowable 
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stress  in  a  bar  element  is  related  to  strain  and  hence  the  relative  displace¬ 


ments  at  each  end  of  the  bar. 

While  this  approach  is  undoubtedly  entirely  valid,  its  extension  to 

other  types  of  finite  elements  may  present  a  major  problem.  Also,  it  is 

believed  that  it  may  require  more  computational  effort  than  the  current 

force  method  approach. 

b.  Force  Method  Background 

With  the  generation  of  optimality  criteria  for  both  stress  and 
displacement  constraints,  the  next  problem  is  the  solution  of  the  resulting 
set  of  nonlinear  equations.  The  most  appropriate  approach  is  through  the 
use  of  some  form  of  Newton-Raphson  technique. 

Attempts  to  solve  the  linearized  form  of  the  total  set  of  governing 
equations  was  relatively  unsuccessful  -  principally  due  to  the  inability  of 
the  search  technique  to  identify  and  distinguish  between  active  and  inactive 
constraints. 

This  difficulty  was  largely  overcome  by  the  introduction  of  a  linear 
programming  technique,  which  provided  an  estimate  of  the  optimal  population 
of  active  constraints.  This  approach  which  is  based  upon  a  local  linearization 
of  the  domain,  leads  to  the  selection  of  full  vertex  corresponding  to  the 
satisfaction  that  the  number  of  constraints  is  equal  to  the  number  of  variables. 

For  stress  constrained  problems  in  which  this  situation  applies  at 
the  optimum,  the  use  of  the  linear  programming  approach  has  been  found 
to  generate  this  optimal  design  in  only  one  or  two  iterative  steps.  The 
validity  of  the  optimality  is  demonstrated  by  checking  that  all  Lagrangian 
multipliers  are  non- negative. 
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Historically,  in  test  problems  involving  only  stress  and  minimum 
size  constraints,  the  optimum  designs  determined  have  appeared  to  be  full 
vertices--i.  e.  ,  the  number  of  constraints  is  equal  to  the  number  of  variables. 
Even  when  it  has  not  been  possible  to  determine  the  optimum  design  rigorously, 
there  have  been  indications  that  the  optimum  design  tends  to  have  a  large 
number  of  active  constraints.  The  question  of  number  of  active  constraints 
is  important  since  it  affects  significantly  the  strategy  required.  Therefore, 
for  a  number  of  variations  on  a  22-bar  truss  theme,  analytic  studies  were 
undertaken  to  determine  the  optimal  design  and  demonstrate  rigorously  the 
active  constraint  populations.  It  was  shown  that,  depending  on  the  values 
of  certain  parameters,  some  optimal  designs  had  22  constraints  and  others 
only  21. 

For  such  problems,  in  which  the  optimum  does  not  occur  at  a  full 
vertex,  the  linear  programming  approach  can  still  play  a  dominant  role. 

The  linear  programming  stage  will  lead  to  a  full  vertex  design  but  one  for 
which  some  Lagrangian  multipliers  are  negative.  These  correspond  to 
constraints  which  should  not  be  active  at  the  optimum  design  and  hence 
should  be  discarded.  For  the  reduced  set  of  constraints^the  Newton- Raphson 
method  can  be  used  satisfactorily.  The  key  to  the  solution  of  the  optimal 
system  is  a  good  reliable  indication  of  which  constraints  are  active  at  the 
optimum  design. 

When  additional  displacement  constraints  are  considered,  the  basic 
approach  is  essentially  unaltered.  Displacement  constraints  appear  even 
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in  the  stress  constraint  case--as  compatibility  conditions.  Thus  the  mathe¬ 
matical  formulation  is  readily  extendable  to  represent  real  external  displace¬ 
ment  constraints.  What  is  different  is  the  uncertainty  regarding  which  of 
these  displacement  constraints  is  active.  By  definition,  all  compatibility 
conditions  are  always  active  constraints;  additional  displacement  constraints 
need  not  be  active.  Also,  displacement  constraints  tend  to  dominate,  some¬ 
times  to  the  exclusion  of  other  types  of  constraints.  Hence,  a  displacement 
constrained  optimum  design  may  have  a  much  sparser  constraint  population 
than  a  full  vertex  design. 

Nonetheless,  the  same  strategy  can  be  applied  as  for  the  purely  stress 
limited  cases.  The  constraint  population  predicted  by  the  linear-programming 
appears  to  be  less  valid  and  greater  reliance  has  to  be  placed  on  constraint 
acquisition  and  discard  algorithms  which  are  incorporated  in  the  full  Newton- 
Raphson  search. 

c.  Report  Organization 

In  Section  2  of  this  report,  a  full  technical  discussion  is  presented 
of  the  various  stages  in  the  development  of  a  force  method  optimization 
program.  Section  2.  1  redefines  the  development  of  optimality  criteria  for 
displacement  constraints  and  highlights  the  problem  of  incorporating  stress 
constraints  in  a  rigorous  manner.  The  fundamentals  of  the  force  method 
are  presented  for  completeness  in  Section  2.2,  since  this  analysis  concept 
is  generally  neglected  in  modern  finite  element  work.  This  discussion  also 
includes  consideration  of  structure  cutters.  Section  2.3  pulls  together  the 
optimality  criteria  approach  and  the  force  method  and  defines  the  overall 


approach  and  governing  equations.  The  problems  and  difficulties  encountered 
in  the  early  work  are  discussed  in  Section  2.4  leading  into  the  development 
of  the  linear  programming  approach  -  Section  2.  5.  The  remainder 
of  the  section  deals  with  the  analytic  study  of  the  22-bar  truss  problem,  the 
concept  of  constraint  discard  and  the  generalized  Newton-Raphson  approach, 
and  finally  the  introduction  of  discrete  displacement  constraints.  Result?  of 
various  optimization  problems  are  inserted  throughout  these  discussions 
where  relevant. 

Section  3  outlines  the  development  and  logic  of  the  computer  program 
which  evolved  from  this  study.  The  final  section  summarizes  the  work  by 
presenting  conclusions  and  recommendations. 
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2.  0  TECHNICAL  DISCUSSION 

2.  1  Optimality  Criteria 

The  concept  of  developing  optimality  criteria  for  structural 
optimization  problems  using  a  Lagrangian  approach  is  well  established. 

This  was  the  basis  of  the  envelope  method  used  in  the  development  of  the 

(2. 

OPTIM  series  of  programs  *  .  A  number  of  variations  on  the  optimality 

criteria  solution  algorithm  have  been  developed,  but  a  comprehensive 
study^  5  ^  has  effectively  demonstrated  that  all  these  algorithms  are 
essentially  similar.  It  must  be  noted,  again,  that  all  these  approaches 
have  been  developed  using  a  vigorous  basis  for  the  consideration  of 
displacement  (and  minimum  size)  constraints.  Stress  constraints,  where 
they  have  been  produced,  have  been  treated  in  a  more  approximate 
manner 

To  understand  this,  it  is  appropriate  to  recapitulate  the 
basic  tenets  of  the  Lagrangian  approach  to  the  development  of  optimality 
criteria. 

The  general  problem  is  the  minimization  of  the  weight  W 
of  a  structure  subject  to  a  number  of  inequality  constraints  of  the  form 
Cj  -  Cj  £  0  acting  on  various  behavioral  response  characteristics  of  the 
structure. 

While  many  types  of  responses  can  be  constrained,  for  the 
purposes  of  this  study  attention  is  confined  solely  to  displacements,  stresses, 
and  minimum  member  sizes.  Other  types  of  constraints  --  stability. 
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dynamic  response,  flutter,  etc.,  are  of  prime  importance  in  structural 
design,  but  they  are  generally  and  genericaliy  related  to  the  three  basic 
constraints  mentioned  above. 

The  standard  approach  is  to  form  a  Lagrangian 

L  =  W  +  Xj  (Cj  -  C.)  (1) 

where  W  is  the  merit  function  (weight)  of  the  system  to  be  minimized. 


C. 

J 


is  the  value  of  some  response  characteristic  which  is 


constrained  to  be  less  than  or  equal  to  an  allowable 


value  Cj 

and  is  an  undetermined  Lagrangian  multiplier. 

Differentiation  of  Equation  1  with  respect  to  the  primary  design  variables 
will  yield  a  set  of  optimality  criteria  which  must  be  satisfied  by  the  optimum 
design.  Differentiation  of  Equation  1  with  respect  to  the  secondary 
variables  (Lagrangian  multipliers)  yields  explicitly  the  constraint  conditions. 
One  major  problem  in  this  or  any  other  approach  to  structural  optimization 
is  the  determination  of  which  constraints  are  active  (e.g.  ,  satisfied  as 
equalities)  and  which  are  inactive  (inequalities).  This  problem  will  be 
addressed  in  a  later  section. 


Under  the  assumption  that  the  constraint  population  can  be 
determined  or  defined,  solution  of  the  combined  set  of  optimality  and 
constraint  equations  (usually  nonlinear)  will  indeed  yield  the  optimum  design. 

To  have  a  Lagrangian,  Equation  1  ,  which  is  capable  of 
explicit  differentiation  for  solution,  it  is  necessary  to  introduce  certain 
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assumptions  and  limitations  on  the  mathematical  model  of  the  system. 

For  the  merit  function,  two  conditions  must  be  satisfied: 

(1)  The  total  weight  is  the  sum  of  the  component  member 
weights. 

(2)  The  weight  of  each  member  is  a  linear  function  of  its 
single  design  variable. 

These  are  fairly  traditional  assumptions  for  weight  minimization 
problems.  The  effect  of  attachments  and  joints  is  ignored  and  attention  is 
restricted  to  essentially  membrane  type  of  behavior.  In  a  finite  element 
idealization,  only  membrane  plates  and  pin-ended  frames  are  used  to  model 
the  structure*. 

The  weight  of  the  structure  can  then  be  written  in  the  form 
m 

w  =  2  W;  A;  (2) 

i  =1 

For  the  constraint  conditions, basically  similar  conditions 
must  be  satisfied: 

(3)  A  constraint  is  applied  to  a  behavioral  characteristic 
of  each  individual  member  of  the  structure  or  on  the 
linear  sum  of  contributions  arising  from  each  member. 

*There  is  an  exception  to  this  wherein  thin-walled  circular  tubes  can  be 
used  to  model  space  frames  with  flexural  members. 
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(4)  Each  constraint  or  constraint  contribution  must  be 


expressible  as  an  explicit  differential  function  of 
the  design  variables. 

Condition  3  clearly  permits  the  consideration  of  member 
stresses  and  minimum  sizes  directly  and  is  consistent  with  the  classical 
virtual  force  method  for  calculating  displacements. 

Condition  4  is  of  much  more  complex  significance  and  has 
a  dominant  influence  on  the  range  of  applicability  of  the  Lagrangian  approach 
to  structural  optimization. 

For  minimum  sizes.  Condition  4  imposes  no  special 
restrictions.  To  be  consistent  with  other  constraints,  it  is  found  to  be 
convenient  to  express  minimum  size  constraints  in  the  form 


£  *  ■ 


i  =  1  ---  m 


(3) 


where  Aj  is  the  current  value  of  a  member  size 

Aj  is  its  minimum  allowable  value, 
and  m  is  the  number  of  elements  in  the  structural  model 

The  corresponding  term  in  the  Lagrangian  is  written 
m 


2  Mf  - ') 

i  =1  1 


where 


/“ 


are  the  associated  Lagrangian  multipliers. 


Displacement  of  a  structure  can  be  computed  using  the  virtual  force  method. 
For  a  pin  jointed  framework,  this  has  the  form 

m 


sf  Li 

EiAi 


(4) 


where 


Sj  are  the  actual  element  forces  due  to  the 
applied  loading 

S'?  are  virtual  element  forces  associated  with  a 


virtual  unit  force  applied  at  the  point  of  the 
unknown  displacement 

Lj  and  Ei  are  the  length  and  moduli  of  the  bar  elements 
For  plate  elements,  similar  expressions  can  be  derived.  The  more  general 
form  of  Equation  4  can  be  written 


<2  =  Spt/  SJ  (5) 

P  i 

where  S  ,  SJ  are  the  vectors  of  element  forces 

and  ^  is  the  assembled  diagonal  element  flexibility 

matrix 

Although  the  displacement  of  Equation  4  has  been  derived  here  through 
what  is  basically  a  force  method  formulation,  this  is  not  an  essential 
requirement.  It  can  be  obtained  with  the  same  degree  of  validity  using 


a  displacement  method. 


In  Equation  4  s|  need  only  be  in  static  equilibrium  with 
the  virtual  load  and  is  hence  independent  of  the  member  sizes  A^.  On  the 

D 

other  hand,  the  basic  force  distribution  is  independent  of  A-'s  only  in 

statically  determinate  structures.  For  more  generally  redundant  systems, 

S^P  is  an  almost  transcendental  function  of  all  member  variables.  As  such, 
the  displacement  constraint  could  not  meet  the  requirement  of  condition  4 
for  differentiation.  In  early  optimality  criteria  studies  it  was  assumed  that 
SiP  was  independent  of  A^  to  overcome  this  obstacle.  The  argument  used 
was  that  the  internal  force  distribution,  even  in  a  redundant  structure, 
was  only  weakly  influenced  by  the  elastic  characteristics  of  the  members 
and  had  negligible  variation  for  small  changes  in  member  size.  Unfortunately, 
this  is  not  generally  true  and,  as  will  be  discussed  later,  the  error  in  this 
assumption  has  hitherto  proved  to  be  the  major  stumbling  block  in  introduction 
of  stress  constraints.  Fortunately,  for  the  case  of  displacement  constraints, 
the  assumption  is  found  to  introduce  no  error.  The  reason  is  that  the 

p 

assumption  of  invariant  ignores  terms  in  the  differentiation  of  the 

type  -  f  SO  .  It  has  been  shown^  5  ^  that  these  terms  do  vanish 

over  the  entire  structure,  since  they  represent  the  virtual  work  of  self- 
equilibrating  systems  and  are  hence  zero. 

Hence  for  displacement,  the  constraint  condition  can  be 
written  in  the  form 


'  <Zli  /  -=■ 

Z_  at  *==  cj 


j  =  1  ---  n 
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where 


ei.  .  "iVh 


are  now  considered  to  be 


‘J 


and 


independent  of  the  variables  A^ 

.th 


Cj  is  the  allowable  value  of  the  j  displacement 
constraint 

The  corresponding  terms  in  the  Lagrangian  are 


n  /  m  x 

Z  X,  drL  ■) 

j=l  i=l  1 


where  are  the  associated  Lagrangian  multipliers. 

For  stress  constraints,  the  form  is  deceptively  simple. 


oi  - 


_£i 

A, 


(?) 


As  discussed  above,  this  will  not  meet  the  requirements  of 
condition  4  directly  for  redundant  systems,  since  Sj  is  a  too  complex 
function  of  all  A^.  There  is  no  corresponding  term  for  stress  constraints 

A  <? 

to  that  which  allowed  the  derivative  —  J  to  be  ignored  in  the  displace - 

dA 

ment  constraint  derivation.  If  it  is  assumed  that  is  invariant,  which 

is  exact  for  statically  determinate  structures  only,  the  use  of  the  optimality 

criteria  approach  will  lead  directly  to  a  fully  stressed  design  --  which  is, 

again,  exact  for  statically  determinate  structures  only. 

The  only  solution  to  this  problem  is  to  include  either  directly 

or  indirectly  terms  which  will  allow  representation  of  L  ^  ^  \  .  One 

\  A  / 
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method  is  to  compute  finite  difference  approximations  but  this  is  generally 
too  expensive  computationally.  Another  approach  lies  in  the  use  of  the 
force  method  as  will  be  discussed  in  the  next  section. 

Omitting  stress  constraints  directly  reduces  the  problem 
to  minimum  size  and  displacement  constraints  only. 

The  Lagrangian  can  now  be  written  as 
m  m  m 


The  third  term  of  Equation  8  for  stress  constraints  can  only  be  strictly  used 
if  the  structure  is  statically  determinate.  The  term  is  included  here  for 
completeness.  Setting  the  differentials  of  Equation  8  with  respect  to 
Ai  equal  to  zero,  leads  to  m  non-linear  equations  with  the  design  variables 
Ai  and  the  Lagrangian  multipliers  and  Xj  as  unknowns.  Selection 

of  a  set  of  active  constraint  conditions  will  define  (a)  which  Lagrangian 
multipliers  are  zero-corresponding  to  inactive  (inequality)  constraints 
and  (b)  provide  a  set  of  k  (where  k£  m)  non-linear  equations  with  design 
variables  as  unknowns. 


If  the  set  of  active  constraints  is  known  accurately  and 
k  =  m,  the  constraint  equations  can  be  solved  using  a  Newton- Raphson 


approach  to  yield  the  optimum  design.  As  will  be  discussed  later  this 
set  of  circumstances  does  not  occur  too  frequently.  Either  k  ^m 
or,  even  more  likely,  the  active  set  of  constraints  is  unknown. 

If  only  one  displacement  constraint  is  active,  the  resulting 
problem,  although  still  nominally  non-linear,  can  be  solved  algebraically. 

This  was  the  basis  for  the  envelope  method ^  ^  from  which  the  OPTIM 

programs  have  been  developed.  In  the  envelope  approach,  each  displace¬ 
ment  constraint  was  considered  in  isolation  and  an  optimum  design  generated 
to  satisfy  that  constraint  alone.  The  interaction  of  several  simultaneously 
active  constraints  was  approximated  by  selecting  a  composite  design  from 
the  envelope  of  the  individual  designs.  This  composite  design  was  modified 
using  a  concept  of  active  and  passive  member  partitioning  wherein  the  design 
of  members  was  allocated  to  be  controlled  by  a  specific  constraint  condition. 
This  partitioning  was  further  modified  by  a  representation  of  member  minimum 
sizes  and  stresses.  The  stresses  were  effectively  represented  by  fully 
stressed  design  criteria.  While  the  envelope  method  has  achieved  some 
success,  its  greatest  weakness  still  lies  in  the  area  of  stress  constraints. 

Alternate  approaches  to  the  approximate  solution  of  the  non¬ 
linear  equations  arising  from  the  differentiation  of  Equation  8  have  been 
very  thoroughly  investigated.  Reference  5  provides  a  comprehensive 
comparison  of  these  methods  of  iterative  solution  of  the  non-linear  problem 
and  attempts  to  identify  the  parameters  which  affect  the  erratic  convergence 
characteristics  frequently  encountered. 
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The  goal  of  the  present  work  was  not  to  duplicate  these 
efforts  but  to  pursue  the  more  elusive  goal  of  incorporating  stress 
constraints  in  a  vigorous  manner.  The  key  to  this  was  felt  to  lie  in  the 
philosophy  of  the  force  method  of  analysis. 

2.  2  Fundamentals  of  Force  Method 

Prior  work  had  indicated  the  use  of  the  force 
method  philosophy  did  provide  the  basis  for  the  vigorous  representation  of 
stress  constraints  in  a  Lagrangian  formulation. 

In  the  vast  majority  of  early  work  in  the  field  of  structural 
optimization,  the  principal  effort  was  directed  towards  the  development 
of  methods  of  redesign  per  se.  The  use  of  a  finite  element  method  of 
analysis  was  considered  merely  as  an  adjunct  for  determining  stresses  and 
displacements  rapidly  and  with  a  minimum  of  computational  effort.  It  was 
stated  that  any  other  method  of  analysis  (e.g.  ,  finite  difference)  would  be 
equally  suitable.  There  was  no  attempt  (or  indeed  any  real  reason)  to 
integrate  the  analysis  and  redesign  philosophies  in  any  real  way.  With 
the  dominance  of  the  displacement  method  of  finite  element  analysis  in  the 
60's  and  70's  this  method  was  fully  developed  and  readily  available  for 
incorporation  into  any  optimization  program  with  a  minimum  of  effort. 

The  force  method  which  had  achieved  initial  prominance  in  the  50' s  was 
effectively  discarded  because  of  some  inherent  disadvantages  it  was  felt 
to  display  in  comparison  with  the  apparently  simpler  displacement  method. 
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In  prior  work  on  optimization^  ^  \  some  use  was  made  of 
the  force  method,  but  this  was  primarily  directed  at  the  reduction  of 
computational  effort  associated  with  iterative  analysis.  The  analysis 
was  not  made  an  integral  part  of  the  otpimization  procedure.  Advantage 
was  taken  of  the  form  of  the  equations  to  be  solved  in  the  iterative  analysis 
to  reduce  computational  times  considerably  compared  with  displacement 
method  techniques. 

Since  there  is  a  general  unfamiliarity  with  the  force  method 
it  is  appropriate  to  discuss  and  define  the  concepts  which  will  be  used  in 
subsequent  sections. 

The  force  method  of  analysis  is  based  upon  the  overall 
enforcement  of  structural  equilibrium  and  the  subsequent  satisfaction  of 
compatibility*. 

For  a  general  finite  element  model  of  a  structure  a  set  of 
overall  equilibrium  equations  can  be  written 

[P]  =  [A]  [S]  (9) 

which  relate  the  externally  applied  loads  [P]  and  the  internal  member 
forces  [S].  If  the  number  of  equilibrium  conditions  is  equal  to  the  number 
of  internal  forces  the  structure  is  statically  determinate  and  [A]  is  a 
square  non-singular  matrix.  If  the  number  of  equilibrium  conditions  is 
less  than  the  number  of  unknown  forces,  the  structure  is  indeterminate 

*In  the  simplest  of  forms  the  displacement  method  can  be  seen  to  be 
the  converse  of  this  concept. 
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or  redundant;  the  degree  of  redundancy  being  defined  by  the  excess  of 
unknown  forces  over  the  number  of  equilibrium  conditions.  It  is  the  latter 
case  which  is  of  significance. 

The  internal  force  distribution  S  is  assumed  to  be  the  sum 
of  two  component  distributions.  The  first  component  is  one  which  is  in 
static  equilibrium  only  with  the  applied  loading  system  while  the  second 
component  arises  from  internally  self-equilibrating  force  systems  of 
undetermined  magnitudes. 

Expressed  mathematically 

S  =  bQ  +  bjX  (10) 

where  bg  is  in  static  equilibrium  with  the  applied  loading 

and  bj  are  values  of  self-equilibrating  unit  force  systems 

whose  magnitudes  X  will  be  determined 
It  is  clear  that  both  bg  and  bj  systems  will  of  necessity,  individually  violate 
compatibility  conditions,  but  their  weighted  linear  combination  will  ensure 
satisfaction  of  compatibility.  The  force  distributions  bo  and  bj  are  not 
generally  unique,  except  in  very  simple  structures.  Their  determination 
may  require  use  of  a  concept  referred  to  as  a  structure  cutter. 

From  the  definitions  of  the  bg  matrix  it  is  clear  that 
there  must  be  zero  forces  in  as  many  members  as  there  are  redundancies. 
These  members  can  be  considered  as  b°ing  cut  in  the  basic  bo  system*. 

#These  cuts  are  not  real  but  are  expressed  as  such  to  give  physical 
illustration  to  the  concept  of  compatibility  violation. 


The  bg  matrix  can  be  written  as 


b  *  P“ 

0  LoJ 


an 


where  bg  are  the  values  of  the  element  forces  due  to  unit  values  of  P. 
In  the  'cut'  members  the  only  forces  will  be  those  due  to  the  redundant 
systems.  If  a  redundant  system  is  assigned  a  unit  magnitude  in  a  'cut' 
member,  the  bj  matrix  may  be  written  as 


,  ] 


(12) 


where  the  partitioning  is  as  in  Equation  11  .  Combining  Equations  11 
and  12  into  10  yields 


[S]  = 


bo' 

rpi 

0 

I 

W 

(13) 


This  matrix  is  square  and  invertible  and  its  inverse  may 


be  written 


L? 

I* 

II 

A1  |  A2 

Lx  J 

;o_ ,  t 

[S] 


(14) 


Equation  (14)  corresponds  to  two  equations,  the  second  of 
which  is  trivial.  The  first  is  exactly  Equation  9  which  defines  overall 
system  equilibrium. 

Thus  to  generate  the  bg  and  bj  matrices  it  is  only 


necessary  to  start  with  the  equilibrium  condition,  Equation  9  .  From  the 
rectangular  matrix  [A]  a  non- singular  square  matrix  [Aj]  is  extracted 


using  a  suitable  rank  technique.  The  bg  and  bj  matrices  are  then 
given  by 

r-v'Azi 

and  bj  =  | - j - J  (15) 

The  extraction  of  fAi]  is  not  a  unique  process  and  can  be 
influenced  by  a  number  of  factors.  The  choice  of  the  bg  matrix  is  of 
minor  significance,  but  the  b^  matrix  will  greatly  affect  the  conditioning 
of  equations  to  be  solved  in  satisfaction  of  compatibility.  A  considerable 
effort^®^  has  been  devoted  to  the  development  of  structure  cutters  which 
improve  conditioning  but  their  study  is  outside  the  scope  of  the  present 
discussion. 


To  enforce  compatibility,  the  magnitudes  of  the  redundant 
forces  must  be  selected  to  ensure  that  the  relative  displacements  of  the 
cuts  in  the  structure  be  set  to  zero.  (An  analogy  to  structural  optimization 
wherein  displacements  are  constrained  to  a  prescribed  value  can  be 
observed.  )  To  accomplish  this  the  virtual  force  method  is  used;  the 
virtual  force  distribution  being  given  by  the  bj  matrix. 

In  matrix  notation  the  compatibility  condition  is  given  by 


b^j-’S  =  0 

bit/bo  +  Q>1  f  bl>  x  =  0 

Hence  X  =  -  (b j1  bj  )~l  (b^  «T  bQ) 


(16) 


(17) 
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Equation  17  is  the  classic  expression  for  the  value  of  the  redundancies. 
While  the  apparent  complexity  of  Equation  17  may  indicate  computational 
co3t,  it  must  be  pointed  out  that  the  order  of  the  matrix  (b^j^  bj)  is  only 
equal  to  the  number  of  redundancies  --a  number  usually  considerably 
smaller  than  the  corresponding  number  of  degrees  of  freedom  encountered 
in  the  solution  of  a  displacement  method  approach. 

Having  discussed  the  force  method  concepts,  it  is  now 
appropriate  to  examine  its  relation  to  optimization. 

2.  3  Role  of  Force  Method  in  Structural  Optimization 

In  Section  2.  1,  the  optimality  criteria  approach  to  structural 
optimization  was  discussed.  The  applicability  of  the  Lagrangian  formulation 
for  minimum  member  sizes,  displacement  constraints  and  stress  constraints 
in  determinate  structures  was  presented. 

In  the  force  method  analysis,  the  internal  force  distributions 
bg,  bj  were  determined  purely  from  static  considerations.  Hence,  apart 
from  the  final  satisfaction  of  compatibility  conditions,  the  entire  analysis 
effectively  treats  only  a  statically  determinate  system.  Reference  to 
Equation  14  indicates  very  clearly  that  if  the  redundant  forces  X  are 
regarded  as  part  of  the  applied  force  set,  the  entire  system  is  indeed 
statically  determinate.  Thus  instead  of  considering  a  redundant  structure 
with  an  external  load  system  P  applied,  the  structure  can  equally  be 
viewed  as  a  determinate  system  with  an  external  loading  system  . 

The  only  additional  requirement  is  that  this  structure  will  be  designed  to 


ensure  that  the  displacements  associated  with  the  (unknown)  forces  X 
shall  be  constrained  to  zero  values.  In  an  optimization  sense  stress 
constraints  may  now  be  considered  because  of  the  static  determinancy. 

Clearly  this  transformation  from  redundancy  to  determinancy 
cannot  be  accomplished  merely  by  making  the  above  statements.  Some 
additional  terms  must  be  introduced  into  the  mathematical  formulations. 

What  has  happened  is  that  the  set  of  design  variables  must  be  expanded  to 
include  the  values  of  the  redundancies  as  unknowns. 

Expressions  for  the  internal  forces  must  now  be  written 
as  linear  functions  of  the  additional  variables  X  (Equation  10)  . 

Condition  4  is  still  satisfied  since  the  new  functional  relationships  are 
still  linear.  The  set  of  optimality  criteria  must  be  expanded  by  the 

generation  of  not  only  — but  also  .  ^ h.  conditions. 

<* A 

The  full  set  of  optimality  criteria  and  constraint  conditions 
applicable  to  a  force  method  formulation  is  presented  in  Fig.  1.  In  the  deriva¬ 
tion  of  the  equations  of  Fig.  1,  two  types  of  displacement  constraints  have  been 
considered,  those  associated  with  the  compability  conditions  and  those  with 
externally  defined  displacement  limitations  (Fig.  9).  The  major  difference 
between  the  two  terms  lies  principally  in  the  absence  of  the  unity  for  the 
compatibility  constraint  value.  Thus  the  stress  constraint  problem  has 
been  solved  effectively  by  transformation  into  a  type  of  displacement  constraint. 
It  is  recognized  that  an  alternate  approach  would  have  been  to  express  each 
individual  (bar)  element  stress  constraint  as  a  relative  displacement 
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Figure  1  -  Optimality  Criteria  and  Constraint  Conditions 


(  9  ) 

constraint  at  the  two  ends.  This  idea  has  been  used  ,  but  it  is  felt  to 
have  the  disadvantage  of  requiring  more  computational  effort  for  reasons 
which  will  be  discussed  in  Section  3.  8,  and  also  its  extension  to  other 
than  bar -type  elements  is  unclear.  The  definition  of  a  simple  relative 
displacement  constraint  criterion  will  not  usually  suffice  for  a  multi-node 
plate  element. 

Having  now  the  rigorous  formulation  which  includes  stress, 
displacement^and  minimum  size  constraints,  a  method  of  solving  the  non¬ 
linear  equations  must  be  derived. 

2.  4  Equation  Solution  and  Constraint  Detection 

In  general,  there  are  potentially  many  more  constraints 
than  variables  and  the  actual  number  of  active  constraints  is  a  small  subset 
of  the  total  number  present.  An  inactive  constraint  is  represented  by  a  zero 
value  of  the  corresponding  Lagrangian  multiplier  and  the  term  hence 
vanishes  from  Equation  8  .  An  active  constraint  corresponds  to  a  non-zero 
multiplier.  Although  no  physical  meanings  are  usually  attached  to  the 
values  of  the  Lagrangian  multipliers,  it  is  clear  that  their  values  are  of 
great  significance. 


Al 

^Ai 


Consider  the  first  two  equations  of  Figure  1,  e.  g.  ,  L  and 

Examination  of  the  form  of  L  and  — reveals  they  can  be  written 

oAi 


the  abbreviated  forms 
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+£  compatibility  terms 
where  represents  all  multipliers 

and  represents  the  remaining  terms. 

By  linear  combination,  the  following  relationship  can  be  formed 


(19) 


L  + 


I  ■  2Lw‘a‘  -I> 


(20) 


For  the  optimum  structure  all  terms  in  L,  except  the  first,  vanish.  Hence 


at  the  optimum  L*  =  W*,  where  an  asterisk  indicates  an  optimal  value. 

b  L 


Also  for  the  optimum  structure 
reduces  to 

W*  =  2W* 

w*  =  > 


=  0.  Therefore  Equation  20 


or 


* 

k 


(21) 


Thus  at  the  optimum  the  weight  of  the  structure  is  given  by  the  sum  of  the 
active  constraint  Lagrangian  multipliers.  * 

While  this  summation  does  indicate  a  physical  meaning, 
the  only  major  use  which  has  been  made  of  Equation  21  has  been  as  an 


*It  is  to  be  noted  that  Lagrangian  multipliers  associated  with  compatibility 
constraints  do  not  contribute  to  summation  of  multipliers  which  equals 
the  weight.  This  is  due  to  the  absence  of  the  unity  which  would  survive 
in  the  combination  of  Equations  19  and  20  . 
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aid  in  determining  when  convergence  on  an  optimal  design  has  occurred. 

Since  the  form  of  the  non-linear  equations  to  be  solved  is 
well  defined,  it  was  felt  that  a  Newton-Raphson  method  to  solution  was 
a  most  practical  approach.  In  a  N-R  formulation,  derivatives  of  all  the 
functions  are  required.  Since  the  entire  approach  via  the  force  method 
has  resulted  in  explicit  terms  of  relatively  simple  algebraic  form,  the 
creation  of  the  derivatives  required  for  the  N-R  solution  is  a  relatively 
straightforward  matter.  The  calculation  of  these  derivative  terms  is 
exact,  which  may  assist  in  convergence,  but  more  importantly,  their 
evaluation  does  not  require  finite  difference  or  similar  calculations  which 
would  involve  repeated  structural  analyses. 

For  the  present  set  of  nonlinear  equations  consisting  of 
optimality  criteria  and  constraints,  the  linearized  N-R  equations  are 
given  in  Figure  2.  All  the  derivatives  are  explicitly  expressible 

in  terms  of  the  constituent  matrices.  The  derivatives  are  listed  in 
Figure  3. 

In  simple  test  problems,  where  the  active  constraint 
population  was  known  from  the  start,  the  N-R  procedure  has  been  shown 
to  converge  rapidly  on  the  optimal  solution.  In  a  small  problem  with  three 
structural  elements  and  two  displacement  type  of  constraints,  the  envelope 
method  required  54  iterations  to  achieve  a  converged  solution.  The  N-R 
approach  converged  in  5  iterations  on  a  design  slightly  lighter  (by  0.  3%) 
than  that  obtained  by  the  envelope  method. 
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FIGURE  2  -  LINEARIZED  NEWTON-RAPHSON  EQUATIONS 


FIGURE  3  -  NEWTON- RAPHSON  DERIVATIVES 


With  the  selection  of  the  N-R  as  the  prime  solution  technique, 


the  major  remaining  task  appeared  to  be  the  generation  of  a  method  for 
differentiating  between  active  and  inactive  constraints. 

Since  all  constraints  are  either  satisfied  as  equalities  or 
their  corresponding  Lagrangian  multipliers  vanish,  the  inequalities 
associated  with  all  constraints  can  be  transformed  into  equalities  by 
expressions  including  their  multiplier.  Thus  a  minimum  size  constraint 
would  be  written 


The  resulting  set  of  equations  tends  to  become  very  large,  but  was  indeed 
a  set  consisting  of  equalities  only.  For  a  problem  with  m  elements, 
p  redundancies  and  constraints  on  both  stresses  and  minimum  sizes, 
a  total  of  3m  +  2p  equations  are  required. 

The  first  problem  to  be  considered  was  the  selection  of  suitable 
starting  points  for  the  linearized  search.  During  the  course  of  this  stage 
of  program  development,  a  considerable  number  of  strategies  were  tested, 
with  very  mixed  results.  The  conclusion  drawn  at  that  time  was  that  the 
search  procedure  was  extremely  sensitive  to  the  initially  selected  values  - 
particularly  the  values  chosen  for  the  Lagrangian  multipliers.  There  was 
some  evidence,  in  a  selected  number  of  cases,  that  areas  tended  to  converge 
on  known  optimal  solutions  fairly  rapidly,  while  multipliers  took  many  more 
iterations  to  stabilize.  From  later  evidence,  accumulated  in  future 
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developments,  it  became  clear  that  these  were  somewhat  erroneous 
conclusions.  Since  the  equations  are  linear  in  the  multipliers,  the  actual 
starting  values  selected  for  the  multipliers  is  immaterial.  What  is  of 
dominant  importance  is  their  population--i.  e.  ,  which  are  zero  and  which 
are  non-zero.  Of  lesser  importance  are  the  starting  values  for  the  member 
sizes. 

Since  these  points  were  not  fully  recognized,  various  schemes 
for  selecting  starting  values  were  investigated.  For  example,  since  the 
sum  of  the  multipliers  was  known  to  equal  the  weight,  at  the  optimum,  the 
following  scheme  was  devised.  It  is  assumed  that  the  applied  load  is  carried 
entirely  in  the  basic  b0  structure,  and  its  members  are  sized  accordingly 
for  stress  limits,  the  other  members  are  set  to  minimum  areas.  Now  for 
each  member,  the  individual  weight  is  computed  and  half  this  value  is 
assigned  to  the  respective  minimum  area  and  stress  Langrangian  multipliers. 
This  approach  assumes,  a  priori,  that  every  potential  constraint  is  active 
and  has  its  own  non-zero  multiplier.  As  an  alternate  form  of  this  starting 
point,  again  the  member  weights  are  used  to  compute  multipliers,  but  the 
full  member  weight  is  assigned  to  either  constraint  dependent  upon  which 
is  assumed  dominant. 

These  approaches  and  others  did  achieve  some  modicum  of 
success  on  three-bar  truss  problems  -  particularly  those  discussed  in  proposed 
contractor  tasks  ,  The  Newton-Raphson  procedure  was  found  to 

be  relatively  unstable^and  extensive  move-limits  had  to  be  introduced  to  prevent 


divergence  occurring  in  early  iterations.  Eventually,  a  limit  of  ±30% 
variation  in  a  design  variable  was  found  to  be  a  satisfactory  limit  for 


smaller  problems.  The  logic  of  the  program  was  designed  to  capture 
additional  constraints,  should  violations  occur  and  to  release  them,  if 
the  associated  multiplier  should  become  negative,  or  if  the  potential 
move  would  take  the  design  away  from  a  constraint. 

The  10-bar  truss  problem  with  the  high  allowable  stress 
in  member  10,  was  next  investigated.  This  proved  to  be  quite  troublesome 
and  required  a  considerable  amount  of  numerical  experimentation  to  obtain 
solutions.  Putting  in  a  design  close  to  the  optimal  did  permit  convergence 
fairly  rapidly,  thereby  demonstrating  overall  correctness  of  the  coding. 

For  cases  in  which  O^q  <.  37^300  psi,  the  optimal  design 
is  fully  stressed  and  the  procedure  did  finally  achieve  the  correct  solutions 
in  a  limited  number  of  iterations  -  less  than  would  have  been  required 
using  a  stress-ratio  method.  For  37500  psi  the  convergence  was 

very  poor  and  sometimes  totally  unstable. 

The  known  optimal  solution  was  reviewed  and  two  potential 
problem  areas  were  tentatively  identified.  Firstly,  there  were  members 
present  which  were  both  fully  stressed  and  at  minimum  area  and  secondly,  the 
optimum  design  had,  by  happenstance  11  active  constraints.  The  former 
condition  was  no  real  novelty  since  it  had  been  encountered  in  some  of  the 
3-bar  test  problems  used  previously.  The  latter  condition,  which  arose 
from  the  geometry  and  loading  of  the  structure  was  felt  to  be  potentially 


more  serious. 


Further  examination  of  the  cases  which  failed  to  converge 
did  indicate  that  the  set  of  linearized  equations  being  solved  were  probably 
singular  -  or  very  close  to  it.  Round-off  accumulation  prevents  exact 
singularity  occurring  but  examination  of  the  values  of  determinants  suggest 
the  presence  of  this  problem. 

The  possibility  of  over-constraint  was  felt  to  be  a  reality  and 
hence  attempts  were  made  to  provide  for  detection  and  elimination  of 
singularities  arising  from  this  source. 

Recourse  was  made  to  the  use  of  the  Cholesky  method  of 
equation  solution  since  the  array  of  equations  is  symmetric.  In  the  triangular 
decomposition  stage,  a  singularity  appears  as  a  zero  (or  near  zero)  term  on 
main  diagonal.  When  such  a  term  is  identified,  the  associated  equation  can 
be  eliminated  from  the  array  by  a  decoupling  process.  This  has  the  effect 
of  winnowing  out  equations  which  are  linearly  dependent  on  any  combination 
of  previous  equations.  Thus  the  order  of  the  equations  will  affect  which 
equations  are  discarded. 

Some  analysis  of  the  algebraic  form  of  the  equations  for 
three-bar  trusses  did  indicate  circumstances  under  which  singularities 
could  occur  -  although  no  generalization  was  possible.  It  was  finally  concluded 
that  it  was  preferable  to  discard  constraints  associated  with  minimum  sizes 
rather  than  stresses.  This  latter  decision  was  purely  arbitrary  and  was 
based  upon  a  judgment  that  stress  constraints  were  of  greater  importance 


than  minimum  sizes. 


The  Cholesky  equation  solver  was  modified  accordingly  to 
remove  singularities  and  the  ordering  of  the  equations  was  changed  to 
ensure  that  stress  constraints  were  favored.  These  changes  had  no  impact 
on  the  small  three-bar  test  problems  but  did  affect  the  10-bar  problems. 

Some  degree  of  success  was  attained  but  the  convergence 
was  still  found  to  be  highly  dependent  on  the  starting  point  selected.  Some 
further  effort  was  expended  on  modifications  to  step  size  limits  but  without 
much  change  in  the  situation. 

Finally,  the  form  of  the  equations  was  reviewed  again  and 
it  was  determined  that  the  strategies  used  had  created  an  undesirable 
situation.  If  any  member  reached  its  lower  limit,  it  would  never  be  possible 
on  any  subsequent  iteration  for  that  member  to  increase  in  size  and  move 
away  from  the  constraint.  That  is,  once  a  minimum  size  constraint  became 
critical,  it  was  retained.  This  was  clearly  unacceptable  and  required  a 
review  of  the  entire  approach. 

From  consideration  of  the  array  of  governing  equations,  it  was 
recognized  that  the  non-linearity  was  confined  to  the  primary  variables  only, 
the  secondary  (Lagrangian)  variables  were  purely  linear.  Although  this  was 
considered  at  some  length,  no  use  could  be  made  of  this  fact  until  Reference  U 
became  available.  This  work  had  the  effect  of  redirecting  the  entire  effort 
into  a  more  fruitful  avenue  of  approach. 


2.5 


Linear  Programmin 


While  the  N-R  solution  technique  combined  with  the  force 
method  does  provide  the  potential  for  the  vigorous  incorporation  of  stress 
constraints,  it  is  clear  that  the  major  problem  of  constraint  detection  still 
remains . 

From  examination  of  the  form  of  the  equations  being 
considered,  it  can  be  seen  that  the  nonlinearities  arise  from  the  design 
variables.  The  equations  are  purely  linear  in  the  Lagrangian  multipliers. 
While  this  suggests  that  some  form  of  decomposition  into  linear  and  non¬ 
linear  problems  might  be  possible,  no  progress  was  made  in  simplifying 
solution  techniques  until  the  work  presented  in  Reference  11  not 
considered. 

In  Reference  11  the  concept  of  duality  is  introduced  and 
is  used  to  linearize  the  problem.  While  the  duality  approach  is  entirely 
valid,  the  same  linearization  can  be  achieved  in  a  more  direct  manner. 

In  Equation  21  it  was  demonstrated  that  the  optimum  weight 
is  the  sum  of  the  Lagrangian  multipliers  associated  with  all  non-zero  value 
constraints.  That  means  the  optimum  weight  is  explicitly  independent  of 
the  (primal)  design  variables.  Hence  if  design  variables  are  selected 
which  satisfy  the  constraint  conditions  either  as  equalities  or  inequalities, 
in  an  arbitrary  manner,  the  optimality  criteria  are  now  a  non-square  set  of 
(m  +  p)  linear  equations  with  n  unknowns  (n>  m  +  p). 


Solution  of  this  problem  can  now  be  accomplished  using  a  standard 


linear  programming  technique  (LPT),  with  the  merit  condition  (weight)  being 
simply  the  sum  o;  the  individual  multipliers.  In  the  LPT,  all  Lagrangian 
multipliers,  except  those  associated  with  the  compatibility  conditions,  will 
be  restrained  to  have  values  greater  than  or  equal  to  zero. 

The  LPT  will  minimize  the  merit  function  by  selecting  which 
variables  shall  have  non-zero  values  (active  constraints)  and  which  remain 
zero  (inactive  constraints).  The  result  is  a  constraint  population,  with  as 
many  constraints  as  there  are  variables. 

With  this  known  constraint  population,  a  new  design  may  be  generated 
readily,  as  will  be  discussed  in  a  later  paragraph. 

This  new  design  may  not,  in  general,  satisfy  all  the  actual  constraints 
in  the  problem  and  hence  may  have  to  be  adjusted,  e.  g.  ,  areas  defined  by 
stress  constraints  may  be  smaller  than  minimum  allowable  sizes.  Dependent 
upon  some  simple  test  criteria,  this  new  design  may  again  be  used  as  the 
starting  point  for  another  iteration  of  the  LPT.  This  cycle  may  be  repeated 
until  convergence  occurs,  in  which  case  the  optimum  design  has  been 
generated. 

It  is  of  interest  to  note  what  will  occur  in  purely  stress  constrained 
problems.  If  other  than  the  force  method  concept  is  used,  all  displacement 
and  compatibility  related  terms  will  vanish  from  the  optimality  criteria 
equations  of  Figure  1.  Examinations  of  the  remaining  terms  in  the  equations 
shows  that  the  array  uncouples  into  m  individual  problems,  each  involving 
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yl A ^  and  only.  The  solution  to  this  then  indicates  that  either  a (minimum 

area)  or  (maximum  stress)  is  present  in  each  member,  i.e.,  the  classic 
FSD  solution.  The  presence  of  the  compatibility  related  terms  associated 
with  the  force  method  formulation  ensures  avoidance  of  this  pitfall. 

The  LPT  approach  generates  a  potential  constraint  population 
at  each  iteration.  It  also  generates  values  for  the  non-zero  Lagrangian 
multipliers.  Experience  indicates  that  these  values  are  of  minor  worth, 
but  the  population  so  defined  is  extremely  relevant.  The  force  method 
approach  is  ideally  suited  to  the  determination  of  a  structure  which  satisfies 
a  prescribed  set  of  constraints.  The  proviso  must  be  made  that  the  number 
of  constraints  equals  the  number  of  variables,  i.e.,  that  the  design  being 
sought  does  lie  at  a  full  vertex  in  the  design  space.  This  is  consistent  with 
the  design  population  prescribed  by  the  LPT  solution  which  will 
have  selected  a  full  vertex  as  optimal  in  the  linearized  design  space.  The 
consequences  of  the  optimal  structure  not  occurring  at  a  full  vertex  will  be 
discussed  in  subsequent  sections. 

In  determining  a  structure  which  corresponds  to  a  given 
(full)  set  of  active  constraints,  the  key  lies  in  the  derivation  of  the  internal 
force  distribution.  An  illustration  of  this  is  the  usual  method  of  generating 
an  FSD  via  stress- ratio. 

An  initial  guess  is  chosen,  the  structure  analyzed,and  areas 
adjusted  to  bring  stresses  to  their  critical  values  and/or  to  their  minimum 
sizes.  The  process  is  simple,  but  many  analysis  iterations  may  be  required 


for  convergence.  This  indeed  is  the  method  conventionally  used  in  conjunction 
with  a  displacement  analysis  technique. 


Assume  now  that  the  distribution  of  critical  constraints  is 
known  a  priori,  i.e.,  which  members  are  fully  stressed  and  which  have 
minimum  areas.  Using  the  force  method  approach,  the  compatibility 
equations  are 


(b^  f  bQ)  +  (b/  f  bt)X  =  0 


’l  A 


Using  a  slightly  different  notation.  Equation  23  can  be  written  in  the  form 


M  i  i  r  ^  ilc 

>1  <V+Ibl  xk>  =  0,  j  =  1  —  n 

1  k 


For  members  in  which  the  stress  is  critical 


bl  Xk 


is  replaced 


by  For  the  other  members  A  is  set  to  Amin.  All  A-values  are  now 

eliminated  and  the  only  unknowns  are  the  redundant  forces  X^.  Solution  of 
the  linear  equations  leads  to  the  required  internal  force  distribution  which 
satisfies  the  prescribed  FSD  conditions --without  iteration.  The  areas  are 
obtained  from  the  known  stresses  and  the  computed  internal  forces. 

If  the  constraint  distribution  given  by  the  L>PT  routine  is 
other  than  FSD-type,  the  equations  to  be  solved  become  nonlinear.  A 
member  may  be  both  fully  stressed  and  have  a  minimum  area.  Consequently 
another  member  size  will  be  completely  unrestrained.  The  presence  of 


active  displacement  constraints  will  similarly  leave  some  members  undefined. 


•4 


iMJtr.  ■ 


In  such  cases,  the  set  of  compatibility  equations  must  be 
expended  by  force-stress -minimum  area  relationships  and/or  displacement 
relationships.  In  either  case  the  number  of  equations  and  variables 
increases  to  (redundancies  +  undefined  areas)  and  the  equations  become  non¬ 
linear.  For  their  solution  a  Newton- Raphson  technique  has  always  proved 
to  be  entirely  adequate  for  problems  in  which  the  defined  constraint  population 
corresponds  to  a  feasible  design.  It  can  occur  that  the  LPT  will  propose 
an  unfeasible  constraint  population.  The  method  of  handling  this  situation 
is  discussed  in  the  program  development  section. 

The  LPT  approach  appeared  to  have  great  merit  for  over¬ 
coming  the  constraint  definition  problems  described  in  the  previous  section. 

A  new  computer  program  was  generated  based  upon  the  LPT 
and  checked  out  on  a  variety  of  small  3-  and  4-bar  truss  problems.  The 
first  major  test  of  the  program's  effectiveness  was  the  classic  10-bar 
truss  problem,  in  which  one  member  has  a  significantly  higher  allowable 
stress  than  the  remaining  nine  members.  This  problem  was  extensively 
discussed  in  Reference  12  ,  for  a  range  of  allowable  stresses  and  had 

been  subsequently  used  as  an  evaluation  tool  for  different  approaches  to 
strength  optimization.  The  new  LPT  approach  was  used  on  this  problem 
and  was  found  to  be  startlingly  successful. 
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2.5.1  Ten-bar  Truss  Problem 


The  ten-bar  truss  shown  in  Figure  4  is  constructed 
of  aluminum  alloy  throughout  (E  =  10  psi,  p  =  0.1  lb/in  ).  The  allowable 
stress  in  all  bars  except  number  10  is  ±25000  psi.  The  allowable  for 
member  10  has  been  varied  in  past  problems  from  ±25000  up  to  ±75000  psi. 

It  has  been  shown^  ^  ^  that  for  —  37500  psi  a  FSD  is  optimum. 

For  C  to  ^  37500  psi  the  optimum  design  (weighing  1497.6  lb)  is  independent 
of  an<3  member  10  is  neither  critically  stressed  or  at  minimum  size. 

Use  of  a  stress  ratio  approach  leads  to  a  totally  different  design  weighing 
1725.24  lb,  an  error  of  approximately  15%  with  a  completely  different 
distribution  of  element  areas. 

The  output  from  the  pilot  computer  program  applied  to  the 
10-bar  problem  with  0  jq  =  ±75000  psi,  is  shown  in  Table  1  .  Computer 
program  development  will  be  discussed  in  a  later  section,  but  it  is  appropriate 
to  mention  here  some  of  the  salient  features  of  these  output  results. 

No  starting  point  design  was  specified  so  the  program 
automatically  selects  minimum  areas  as  the  default  option.  Experience 
has  indicated  that  such  designs  may  be  very  poor  approximations  of  optimal 
systems,  hence  one  stress-ratio  cycle  is  optionally  permitted  to  generate  a 
more-or-less  feasible  design.  This  design  is  still  very  different  from  either 
the  (known)  optimum  or  the  converged  FSD.  This  cycled  design  is  used  in 
the  LPT  routine  and  results  in  the  proposed  population  of  minimum  areas  and 
fully  stressed  members  shown  in  Table  1.  The  negative  sign  on  stressed 
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/INPUT 

/JOB  TIME*! ,60 
/ INCLUDE  RGHSC 
/INCLUDE  NEWT 
/INCLUOE  SU8S 
/INCLUDE  BARTEN 

/mm 

PAUSE  ENTER  T0L.T0L1 , IRO.ISR.IPRINT.KOOEI ,K0DE2 
I .£'8,1 .E-7,0. 0,-3 -0,0 


STRESS  RATIO  AREAS 


I  0.81854E  01 

2 

0.23950E  01 

3  0.16050E  01  4 

0.16050E  01 

5 

0.78146E  01 

6  0.14196E  01 

7 

0.53946E  01 

8  0.59190E  01  9 

0.22698E  01 

10 

0.11290E  01 

X 

I  -0.2J033E  05 

2 

0.30183E  05 

MEM  1234 

5 

6  7  8  9  10 

AREA  0  0  Q  I 

0 

1  0  0  0  0 

STR  -1-1  1  1 

1 

0-1  1-1  0 

X  VALUES  1  -0. 

18533E  05  2 

0.47683E  05  3 

0.12D12E  01 

X  VALUES  1  -0. 

12695E-0!  2 

0.97656E-03  3 

0.17090E  00 

N-R  AREAS 

1  0.81000E  01 

2 

Q.3900QE  01 

3  0. 10000E  00  4 

0.10000E  00 

5 

0.79000E  01 

6  0.10000E  00 

7 

0.55154E  01 

8  0.57983E  01  9 

0.14142E  00 

10 

0.36769E  01 

X  VALUES  1  -o.; 

25000E  04  2 

0.25000E  04 

LAGRANGIAN  MULTIPLIERS 

1  0.0 

2 

0.0 

3  0,0  4 

0. 60001 E  01 

5 

0.0 

6  0.36000E  01 

7 

0.0 

8  0.0  9 

0.0 

10 

0.0 

11  0.33400E  03 

12 

0.78000E  02 

13  0.66000E  02  14 

0.60000E  02 

15 

0.22200E  03 

16  0.0 

17 

0.15600E.03 

18  0.42000E  03  19 

0.13200E  03 

20 

0.0 

21  0.69333E  02 

22 

-0.69333E  02 

ITER  NO.  1 

SUM  MEM  WGT 

» 

0.14976E  04 

SUM  LAM  •  0.14976E  04 

M.0070  END 

M.0072  BEGIN 

TABLE  1  - 

10-8AR  TRUSS  -  COMPUTER  RESULTS 
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members  indicates  a  compressive  allowable.  It  is  to  be  noted  that  member 
10  is  undefined,  while  member  3  is  both  fully  stressed  and  of  minimum  area. 

A Newton-Raphson  routine  is  used  to  generate  a  design 
corresponding  to  the  specified  constraints  using  the  previous  design  as  a 
starting  point.  Since  the  structure  is  doubly  redundant  and  one  member 
area  is  undefined,  there  are  only  three  nonlinear  equations  and  three  variable 
(Xj,  X2  and  Ajq).  Only  three  iterations  are  required  for  convergence. 

With  these  new  values  of  the  redundant  forces,  the  known  stresses,  minimum' 
sized  members  and  an  area  for  member  10,  the  new  design  is  generated  and 
printed  out.  It  can  be  seen  that  member  4  is  at  minimum  size  but  this  is 
not  recognized  as  an  active  constraint.  This  is  a  rare,  but  not  unique, 
situation  which  arises  from  the  geometry  and  loading  at  Gridpoint  3. 

Members  3  and  4  are  effectively  linked  together.  The  presence  of  a  surplus, 
unused  constraint  has  no  influence  on  the  problem  or  solution  technique. 

A  shortage  of  constraints,  on  the  other  hand,  is  of  major  consequence  as 
will  be  discussed  later.  For  this  design,  the  corresponding  non-zero 
Lagrangian  multipliers  can  be  obtained  from  the  linear  solution  of  the 
optimality  criteria  equations.  Of  the  printed  set,  the  first  ten  multipliers 
are  those  associated  with  minimum  sizes,  the  next  ten  with  stresses  and 
the  final  two  are  the  compatibility  constraints. 

Since  all  multipliers  associated  with  non-zero  type  constraints 
have  positive  values,  the  design  is  optimal.  As  an  additional  check  the  sum 
of  the  multipliers  does  equal  the  weight.  This  solution,  with  its  extremely 
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rapid  convergence  is  an  excellent  demonstration  of  the  power  of  the  L.  P.  / 
Force  method.  The  number  of  steps  was  minimal  as  was  the  computational 
effort. 

The  next  stage  in  the  check-out  process  was  the  application 
of  the  program  to  a  larger  problem  --  the  22-bar  double  truss. 

2.  6  Twenty-two  Bar  Truss  Problem 

The  22-bar  truss  problem  is  an  extension  of  the  10-bar 

problem  insofar  as  it  consists  of  two  10-bar  cantilever  frames  connected 

at  their  extremities  by  a  pair  of  hangars  (Figure  5  ).  A  single  load 

of  100,000  lbs  is  applied  vertically  at  the  extremity  of  the  hanger.  While 

7 

the  modulus  is  held  constant  at  10  psi  throughout,  various  allowable  stresses 
and  material  densities  were  to  be  considered.  The  truss  problem  has  proved 
to  be  a  significant  one  and  has  been  investigated  by  a  number  of  researchers^’ 
Due  to  a  desire  to  explore  the  effects  of  varying  critical  parameters  and  also 
as  a  result  of  confusion  as  to  specified  values  of  some  parameters,  a  number 
of  cases  have  come  into  being.  This  confusion  has  been  compounded,  possibly 
by  the  fact  that  all  problems  have  been  labeled  either  Case  A  or  Case  B. 
Whatever  particular  values  for  stress  allowables,  densities  or  minimum 
sizes,  the  essence  of  the  problem  here  is  to  demonstrate  the  manner  in  which 
the  applied  load  is  distributed  optimally  as  the  strength/density  ratios  are 
varied  between  the  two  half  (10-bar)  frames.  The  stress -ratio  method  is 
clearly  useless  for  this  type  of  problem  since  it  includes  no  reference  what¬ 
ever  to  material  density- -clearly  a  major  factor. 


Some  attempts  were  made  to  run  some  cases  using  the 
program  developed  previously.  The  results  obtained  were  extremely 
variable.  Some  cases  produced  solutions  readily,  whereas  others  failed  to 
converge  at  all.  This  led  to  a  re-evaluation  of  the  concepts  used.  In 
particular,  the  question  of  the  number  of  active  constraints  at  the  optimum 
design  was  considered.  In  all  stress  constraint  test  problems  which  could 
be  reviewed  either  the  number  of  active  constraints  was  equal  to  the  number 
of  variables  (full  vertex)  or  it  was  impossible  to  state  reliably  how  many 
constraints  were  active.  This  latter  situation  occurs  in  solutions  where, 
for  example,  many  stresses  or  areas  are  within  a  few  percent  of  their 
critical  values  and  there  is  no  overt  determination  of  the  active  constraint 
population. 

Since  this  question  of  number  of  active  constraints  is  of 
prime  consideration,  a  more  detailed  study  of  the  22-bar  problem  was 
deemed  appropriate.  Although  this  structure  (in  its  various  forms)  has 
been  difficult  to  optimize  using  general  purpose  optimization  methods,  the 
generation  of  optimal  solutions  analytically  is  relatively  straight  forward. 

In  all  the  cases  of  interest,  the  allowable  stresses  and 
densities  are  uniform  throughout  each  10-bar  truss  which  form  the  two 
principal  halves  of  the  structure.  The  approach  to  the  analytical  optimization 
is  through  substructuring.  The  individual  10-bar  trusses  are  optimized 
parametrically  in  terms  of  the  applied  load  and  the  two  halves  are  coupled 
through  a  compatibility  condition.  This  approach  is  only  possible  here  due 


to  the  special  geometry  of  the  structure,  although  the  substructuring  concept 

(15 

has  been  proposed  for  general  use  in  the  optimization  of  complex  structures 
From  the  work  performed  in  the  study  of  the  10-bar  truss, 

discussed  previously,  it  has  been  conclusively  demonstrated  that  a  fully 

constrained  design  is  optimal  for  uniform  material  and  strength  properties. 

For  a  vertical  load  applied  at  the  outermost  lower  vertex,  the  optimal 

population  has  been  determined.  The  population  is  that  which  normal 

engineering  intuition  would  suggest,  inasmuch  as  there  is  little  load  diffusion 

from  principal  determinate  load  paths  when  the  redundant  elements  are 

introduced. 

To  optimize  the  structure,  the  force  method  approach  is  used. 


a)  Basic  Structure 


b)  Determinate  (bg)  System 


c)  Redundant  (bj)  Systems 


FIGURE  6 


10-BAR  TRUSS 


The  basic  structure  is  shown  in  Figure  6  a).  The  constituent  determinate 

bg  and  redundant  bj  systems  are  shown  in  b)  and  c).  The  choice  of  these 
systems  is  non-unique,  but  the  members  retained  in  the  bQ-system  are  those 
which  form  the  natural  primary  load  bearing  structure.  Provided  the  ratio 
between  the  applied  load  and  the  stress  allowable  is  sufficiently  high,  the 
members  1,  3,  4,  8  and  9  present  in  the  bo-system  are  fully  stressed. 

The  remaining  members  2,  5,  6,  7  and  10  have  minimum  areas.  It  is 
possible  that  the  applied  load  could  be  so  small  that  some  or  all  of  the  primary 
bQ-system  members  would  have  to  be  set  to  their  minimum  areas  instead  of 
being  fully  stressed.  In  the  present  analysis,  this  situation  has  not  been 
found  to  occur.  For  the  systems  shown,  the  relevant  matrices  are 
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where  f  is  the  diagonal  matrix  of  unit  sized  element  flexibilities 
L  is  the  length  of  a  vertical  or  horizontal  element  (  =  360") 
and  E  is  the  uniform  elastic  modulus. 

With  the  information  that  members  1,  3,  4,  8,  and  9  are  fully  stressed  to 
the  allowable  value  and  the  remaining  members  have  a  minimum  area  A*, 
the  usual  compatibility  conditions  can  be  written  as 
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(26) 


It 

C\J 

X 

+ 

r— 

X 

+ 

tt 

CM 

-  2  <TA*  \ 

/ 

x1  +  ( 3+2JT)  x2  = 

3  <TA*  / 

Hence 

Xt  *  -.68767264 CT"A * 

^cta*  , 

& 

X2  =  .63270460  <TA* 

<*2  <Ta*  1 

The  resultant  internal  force  distribution,  S,  is 


S  = 


2P  +  ^  0“a* 

*2aV 

-P  +  o<]  (T A* 

-p  +  a* 

+  e(2)CTA* 
^2(TA* 

-  J2‘<1  CTa* 

-J7p  -  A* 

J2P  -  J7o/2(Ta* 

-jt^2o-a*  . 


The  correspond!' ng  areas  are 


A  = 


2  p/^  +  o<1  a* 

A* 

p/cr  -  a* 

P/(T-  c<2  A* 

A* 

A* 

A* 


r2(P/ir+^1  a*) 

2(P^--^2  A*) 
A* 
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(29) 


The  weight  of  the  optimized  10-bar  truss  is 


W  *  [8P/cr  +  ( 3+2 12  +  2  -<2)A*3  360  ^  (30> 

where  is  material  density. 


Thus  the  weight  of  each  half-structure  is  expressible  purely  as  a  function 
of  the  load  P  at  the  free  end.  Since  the  total  load  on  the  structure  is  known, 
the  problem  is  now  reduced  to  a  single  variable  optimization  which  defines 
the  sharing  of  the  externally  applied  load  between  the  left  and  right-hand 
halves. 

The  problem  is  further  simplified  when  the  deflection  of  each  half 
is  computed.  Using  the  unit  virtual  load  method,  only  the  members  in  the 
determinate  bQ  system  need  be  considered.  Since  these  members  are  by 
definition  fully  stressed,  their  strains  are  independent  of  their  areas, 
hence  the  total  deflection  is  independent  of  the  applied  load. 


1  .e . 


S 


*  zL0'  *  bi 

1 ,3, 4, 8, 9 


(31) 


Note  here  that  <3^  is  intended  to  signify  both  positive  and  negative  stresses 

which  must  be  accounted  for  in  performing  the  summation. 

Using  the  values  of  Equation  25  the  deflection  of  the  10-bar  frame  is 

S  =  (32) 

E 

Since  L  and  E  are  not  variables  in  this  problem,  the  deflection  of 
each  half-frame  is  directly  proportional  to  the  allowable  stress  only.  Thus 


when  the  allowable  stresses  are  different  on  each  side,  the  higher  stressed 


side  will  always  deflect  more,  irrespective  of  which  side  carries  the 
predominant  share  of  the  applied  load. 

In  any  particular  case,  the  stresses  are  specified  and  the  individual 

deflections  are  known  (e.  g.  ,  and  S  ).  The  ends  of  the  hangars 

J-  R 

(members  21  and  22)  are  joined  together.  Hence  the  difference  in  stress 
between  the  two  hangers  is 

^0^  =  —  (  -  Cf  R)  (3  31 

where  I /  is  the  hanger  length  (  =  3600")  and  the  subscripts  indicate  hanger, 

LHS  and  RHS,  respectively. 

Finally,  from  an  inspection  of  the  problem  parameters  (allowable 
stress  and  density),  it  can  be  readily  decided  which  half  structure  should 
carry  the  major  portion  of  the  applied  load.  A  logical  corollary  from  this 
is  the  fact  that  the  hanger  on  the  lightly  loaded  side  will  have  a  minimum  area. 
It  follows  therefore  that  the  optimum  structure  will  have  at  least  21  active 
constraints.  The  existence  or  otherwise  of  the  22nd  depends  on  the  values 
of  the  parameters. 

To  complete  the  optimization,  it  is  assumed  that  the  LHS  always  has 
the  higher  allowable  stress.  Initially  it  is  taken  that  the  stress/density  ratio 
ensure  that  the  major  load  is  carried  on  the  LHS.  Substituting  numerical 
values  into  Equation  30  the  weights  of  the  two  halves  can  be  written  as 
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WL  +  2 . 554968A*  J  360  /°L 

WR  =>  +  2.554968A*3  360 

5 

The  applied  load  of  the  10  lb  is  the  sum  of  the  individual  hanger 
loads,  i.e. 


(34) 


P  +  Pp  =  io- 
L  K 


Hence 


r  ®pR 

“*■  ■  L-  n 


8x10  *“1 

+  — +  2 . 554968A  360 

or  l  J 


A 


(35) 


(36) 


*  * 

Since  the  major  load  is  on  LHS,  Aj^  =  (Ajj  =  minimum  area  of 

hanger)  and 


H 


'H 


(37) 


Because  the  LHS  deflects  more  than  the  RHS,  the  stress  in  the  RH  hanger 
is  less  than  that  of  the  LH  hanger. 


Then  C"  l 


H 


&  RH 

■Ait'h 

PR 

A* 

<  <4  -  <A> 

L 

area  A 

^  is  given  by 

(  P_l_ 

105  -  Pr 

^LH 

Pr  .  SJL  (Cr 
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The  total  weight  is 


WT0T  *  WL  +  WR  +  WLH  +  WRH  <40> 

While  this  expression  is  somewhat  complicated  algebraically,  substitution 
ot  numerical  values  for  known  parameters  simplifies  it  to  the  form 


“tot  •  C1  ♦  c2  PR  *  . — — —  (41) 

(Pr-c4) 

where  C],  C£,  C 3  &  C4  are  numerical  coefficients 
The  weight  of  the  structure  is  hence  a  nonlinear  function  of  the  single  variable 

Pp^.  The  minimum  weight  may  then  be  either  a  free  minimum  from  Equation 

41  or  may  be  constrained.  Differentiating  Equation  41  WRT  PR  yields 


d WTQT  _  c  _  c3 

a  pr  2  (Pr-C4)2 


Since  the  RH  hanger  has  minimum  area,  if  its  stress  (Equation  37  ,  computed 
using  PR  given  by  Equation  42  ,  is  greater  than  its  allowable,  then  the  RH 
hanger  is  at  both  minimum  area  and  maximum  stress  and  Equation  42  does 
not  apply.  In  this  case 
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and  the  optimum  design  has  22  active  constraints.  With  the  determination 
of  PR  from  either  Equation  42  or  Equation  43  ,  the  optimum  design  is 


1 


completely  defined.  Both  fully-  and  partially-constrained  designs  can 
therefore  exist  for  stress  constrained  problems.  If  the  stress /density 
ratio  transfers  the  loading  to  the  RHS,  a  similar  procedure  is  followed 
to  obtain  the  optimum  designs.  Again,  depending  on  the  values  of  the 
parameters  both  21  and  22  active  constraint  designs  may  occur. 

During  the  study  of  this  basic  problem  by  various  investigators,  a 
number  of  different  combinations  of  parameters  have  been  used,  labeled 
either  Case  A  or  Case  B. 

To  attempt  to  bring  some  clarity  into  the  situation,  the  eight  cases 
given  in  Table  2  are  defined.  It  is  believed  that  these  eight  cases 
encompass  the  principal  problems  studied  by  various  investigators.  It  is 
left  to  the  reader  to  correlate  the  new  and  old  designations. 


TABLE  2 

22  BAR  TRUSS  -  COMPUTER  RESULTS 


CASE 

<r 

ksi 

LHS 

P 

RHS 

P 

A* 

<r  . 

ksi 

HANGER 

P 

A* 

1 

50 

.1 

.001 

25 

.1 

.001 

5000 

.1 

.001 

2 

30 

.3 

.001 

25 

.1 

.001 

5000 

.1 

.001 

3 

50 

.1 

.001 

25 

.1 

.001 

500 

.1 

.001 

4 

30 

.3 

.001 

25 

.1 

.001 

500 

.1 

.001 

5 

50 

.1 

.01 

25 

.1 

.01 

5000 

.1 

.001 

6 

30 

.3 

.01 

25 

.1 

.01 

5000 

.1 

.01 

1  7 

50 

.1 

.01 

25 

.1 

.01 

500 

.1 

.001 

8 

30 

.3 

.01 

25 

.1 

.01 

500 

.1 

.001 
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Using  the  analytical  approach  described  above,  the  following  optimal 
solutions  have  been  determined  (Table  3). 


TABLE  3 

22  BAR  TRUSS  OPTIMAL  SOLUTION  OF  PRINCIPAL  CASES 


— 

CASE 

r . r  - 

2 

3 

4  i  5 

6  7 

i  8 

WEIGHT 

605.096 

1202. 1 74 

654.04 

1232.94 

606.752 

1312.764  655.70 

1236.2 

NO.  OF 

1 

t 

i 

CONSTB. 

21 

21 

22 

22 

_ - 

21 

21  |  22 

|  22 

With  the  determination  that  stress  constrained  problems  could  exhibit  less 
than  fully  constrained  optimal  designs  it  was  necessary  to  revise  the 
previously  developed  approach  to  accommodate  this  situation. 
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Constraint  Discard  and  Generalized  Newton-Raphson 


With  the  determination  that  optimal  stress  limited  designs 
could  be  less  than  fully  constrained,  it  became  necessary  to  revise  the 
LPT  approach. 

From  the  discussion  presented  in  Section  2.4,  it  is  clear 
that  the  operation  of  a  successful  search  procedure  is  highly  dependent 
upon  some  prior  knowledge  of  active  constraint  population.  The  LPT  will, 
by  definition,  always  select  a  full  vertex  in  the  linearized  design  space.  If 
the  LPT  approach  is  valid,  the  optimal  population  should  be  a  subset  of 
the  LPT  population.  Hence,  one  approach  is  to  generate  a  feasible  (but 
non-optimum)  full  vertex  design  and  then  discard  unwanted  constraints  in 
some  rational  manner.  If  a  less  than  fully  constrained  optimum  design  is 
used  as  an  input  guess  to  the  LPT  routine,  a  full  complement  of  constraints 
will  be  identified,  but  the  Lagrangian  multiplier  associated  with  the  spurious 
constraint  will  be  zero.  Unfortunately,  this  does  not  prove  to  be  of  much 
assistance,  since  the  spurious  multiplier  generated  by  the  LPT  will  be 
non-zero  when  an  off-optimum  design  is  used  as  input.  Also  since  experience 
indicates  that  the  values  of  the  multipliers  vary  by  several  orders  of 
magnitude  in  any  one  problem,  it  is  not  practical  to  use. absolute  or  relative 
size  as  a  criterion  for  constraint  discard. 

After  some  experimentation  with  a  variety  of  techniques, 
which  included  arbitrary  or  preferential  discard  of  constraints  associated 
with  minimum  size  members,  a  new  strategy  was  developed. 
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In  the  LPT  routine,  the  variables  (Lagrangian  multipliers) 
are  restrained  to  have  values  greater  than  or  equal  to  zero.  Negative 
multipliers  are  not  permitted  (except  those  associated  with  compatibility 
constraints).  Subsequently  a  new  design  is  created  which  has  active 
constraints  corresponding  to  the  population  generated  by  the  LPT  routine. 

For  this  new  design,  a  new  set  of  Lagrangian  multipliers  can  be  generated 
by  linear  solution  of  the  optimality  criteria  equations. 

If  all  the  multipliers  are  positive,  the  design  is  optimal. 

The  presence  of  one  or  more  negative  multipliers  indicates  that  non-optimal 
constraints  associated  with  these  negative  multipliers  leaves  a  reduced  set 
of  potential  active  constraints  for  the  optimum  design. 

As  discussed  in  Section  2.4,  the  use  of  a  Newton-Raphson 
technique  was  investigated  for  the  solution  of  the  full  set  of  nonlinear 
equations  consisting  of  the  optimality  criteria  and  all  the  potential  constraints. 
This  approach  failed  because  of  the  difficulty  of  identifying  meaningful  sets 
of  active  constraints.  If  an  incomplete  set  was  selected,  there  was  a 
recurrent  problem  with  singularity. 

In  the  present  approach,  the  number  of  active  constraints  is 
specified  a  priori.  Also  the  number  of  constraints  is  usually  a  relatively 
high  proportion  of  the  number  of  variables --thereby  seemingly  avoiding 
the  problems  encountered  previously. 

The  governing  nonlinear  equations  were  developed  in  Section 
2.4,  along  with  the  Newton  linearization. 
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When  the  constraint  population  is  known  the  full  set  of 
Newton  equations  must  be  reduced  to  provide  the  non-singular  set  to  be 
solved. 

If  the  structure  has  m  elements,  each  of  which  has  both 
stress  and  minimum  size  constraints  and  there  are  n  redundancies,  the 
total  number  of  equations  available  is  (3m  +2n).  The  Newton  equations 
can  be  arranged  in  the  symmetric  array  shown  in  Figure  7. 

If  a  member  is  known  to  be  at  a  minimum  size  constraint,  it  is  not  strictly 
a  variable.  Hence  the  corresponding  optimality  criterion  does  not  apply 
and  can  be  deleted  along  with  the  associated  variable,  Deletion  of 

the  variable,  £  A^,  will  result  in  the  elimination  of  the  constraint  equation 
associated  with  that  minimum  area.  The  minimum  area  constraint  equations 
for  unconstrained  members  do  not  apply.  Hence  the  complete  set  of  minimum 
area  constraint  equations  are  deleted.  Finally,  constraint  equations  and 
Lagrangian  multipliers  for  non-stress  constrained  members  are  not  applicable. 
This  reduces  the  number  of  applicable  equations  to  the  more  manageable 
order  (2n  +  m  +  s  -  a)  where  s  is  number  of  stress  constraints  and  a  is 
number  of  minimum  sized  members. 

When  a  feasible  design  has  been  generated  by  the  LPT  section, 
the  corresponding  Lagrangian  multipliers  are  determined.  Constraints 
corresponding  to  any  negative  multipliers  are  discarded  and  the  surviving 
nonlinear  equations  are  solved  iteratively  using  the  full  Newton-Raphson 
method. 
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FIGURE  7.  SYMMETRIC  ARRAY-NEWTON  EQUATIONS 


(2)  For  each  element  without  stress  constraint,  fs.  equation  and  <5v.  variable  are  deleted. 


As  starting  points,  the  current  values  of  member  sizes 
and  redundancies  are  used.  Since  the  equations  are  linear  in  the  Lagrangian 
multipliers,  the  Newton  solution  will  always  generate  the  same  final  values 
for  the  multipliers,  at  each  step,  entirely  independent  of  initially  selected 
guess  values.  The  rate  of  convergence  is  therefore  unaffected  by  selection 
of  the  starting  point  for  the  multipliers.  It  is  controlled  by  the  definition  of 
the  constraint  population  and  is  also  influenced  by  the  initial  guesses  for 
the  areas. 

Also,  as  is  common  in  many  linearized  solutions  of  non¬ 
linear  problems,  the  early  iterations  have  to  be  controlled  to  prevent 
instabilities  occurring.  The  usual  method  used  here  is  the  imposition 
of  move-limits  on  variables.  These  strategies  were  used  in  the  program 
development. 

To  test  the  extended  computer  program,  yet  another  three- 
bar  truss  problem  was  devised.  This  one  is  relatively  unique  in  that  it  has 
only  two  active  constraints.  Details  of  the  problem  are  given  in  Figure  8. 

The  program  had  no  difficulty  in  determining  the  optimal  design. 

With  the  development  of  the  extended  program,  the  22-bar 
structures  which  had  contributed  so  much  to  the  detection  of  the  non-fully 
constrained  problems,  could  not  be  handled. 

Of  the  eight  problems  tested,  all  except  Cases  1  and  3  executed 
successfully.  Table  4  summarizes  the  numbers  of  steps  required  for 
convergence.  In  all  a  single  stress-ratio  step  was  used  initially  to  produce 
a  more-or-less  feasible  design. 


As  can  be  seen  from  Table  4  the  results  generated  correspond 


exactly  to  those  given  in  Table  3. 

Although  no  solution  could  be  generated  for  Cases  1  and  3 
using  the  general  starting  point,  the  program  did  recognize  the  optimal 
designs  when  supplied  as  input.  Table  5  provides  details. 

TABLE  4 

22  BAR  TRUSS  -  ITERATION  HISTORY 


STRESS  NO.  OF 

RATIO  NO.  OF  L-P  NO.  OF  N-R  FINAL  ACTIVE 

ITERATION  ITERATIONS  ITERATIONS  WEIGHT  CONSTRAINTS 


REMARKS 


1202.1  7 

21 

— 

1232.9 

22 

606.75 

21 

1312.76 

21 

655.70 

22 

1236.2 

22 

L-P  Stages 
Di verged 

Final  Solution 
is  linear 

Fastest  Con¬ 
vergence 


Slow  N-R 
Convergence 


Final  Solution 
is  linear 


,  TABLE  5 

22  BAR  TRUSS  -  TEST  CASE  RESULTS  USING  OPTIMAL  INPUT 


CASEfSTRESS  NO.  OF  L-P  NO.  OF  N-R  FINAL  ! NO .  OF  ACTIVE 
NO.  RATIO  ITERATIONS  ITERATIONS  WEIGHT  CONSTRAINTS  REMARKS 


605.094 


654.04 


Opti ma 1 
Sol uti on 


V: 


Before  assessing  the  reasons  for  the  difficulties  with  Cases  1  and  3,  it  is 
appropriate  to  evaluate  some  of  the  computational  requirements  of  the  new 
procedure. 

The  22-bar  truss  has  5  redundancies.  Hence,  each  analysis 
requires  solution  of  5  equations.  The  stress-ratio  requires  one  such 
analysis. 

The  LPT  stage  extracts  the  optimal  solution  from  a  set  of 
(22+5)  =  27  equations  involving  (2x22+5)  =  49  unknowns.  Typically,  the 
solution  requires  40-50  row/column  interchanges.  A  double  precision 
Bell  LPT  routine  is  used.  The  double  precision  was  introduced  because  of 
an  early  impression  of  conditioning  problems.  There  is  no  evidence  that 
double  precision  is  necessary  for  the  cases  actually  solved.  The  reduced 
N-R  used  between  LPT  stages  requires  the  iterative  solution  of  (5+NIQ) 
equations,  where  NIQ  is  the  number  of  doubly  defined  elements.  When  the 
solution  is  linear,  5  equations  are  solved  once. 

In  the  major  N-R,  the  number  of  symmetric  equations  solved 
was  only  31,  since  10  stress  constraints  and  11  minimum  areas  were  defined 
for  the  optimal  population.  The  generation  of  all  the  terms  in  the  array 
solved  is  simple,  direct  and  requires  no  auxiliary  analyses  to  generate 
approximate  derivatives. 

Turning  to  the  two  cases  on  which  the  program  failed,  there 
appears  at  first  glance  to  be  no  obvious  reason  for  the  inability  of  the  LPT 
routines  to  generate  feasible  constraint  distributions.  It  is  true  that  these 


cases  represent  the  most  radical  variations  in  allowable  stresses  coupled 
with  extremely  small  minimum  areas.  Examination  of  the  values  of  the 
Lagrangian  multipliers  does  indicate  a  potential  conditioning  problem. 

The  sum  of  the  Lagrangian  multipliers  is  the  merit  condition  for  the  LPT 
In  Case  1,  the  multipliers  range  in  value  from  140.  31  down  to  .00276  - 
which  immediately  suggests  a  potential  problem.  In  the  successful  Case  2, 
the  corresponding  range  is  283.  77  to  .  00276  -  and  no  difficulty  was 
encountered.  It  is,  nevertheless,  entirely  possible  that  a  conditioning 
problem  does  exist  in  the  LPT  stage.  Switching  to  double  precision  in  the 
LPT  routine  may  be  only  half  the  solution.  It  may  be  necessary,  for 
problems  as  extreme  as  those  considered  herein,  to  formulate  the  terms 
in  the  LPT  in  double  precision. 

Finally,  while  comparisons  are  odious,  it  is  relevant  to 
compare  the  present  solutions  for  22-bar  trusses  with  a  few  available  from 
alternate  sources. 

The  22-bar  was  discussed  extensively  in  Reference  9. 

Two  cases  were  presented.  From  the  results,  these  are  believed  to 
correspond  to  Cases  1  and  2  of  the  present  report.  Case  1  yielded  a  final 
weight  of  606  lb.  and  Case  2  a  weight  of  1203  lb.  Both  required  a  large 
number  of  iterations  to  converge.  For  Case  2,  an  even  lower  weight  of 
1188  lb.  is  quoted,  but  it  is  not  clear  as  to  whether  this  represents  a  truly 


The  team  of  Bartholemew  and  Morris  at  RAE  has  developed 
the  SCICON  structural  optimization  system^^  and  applied  it  to  the  22-bar 
problem^  For  Case  2  a  weight  of  1204.  38  lb.  was  generated  in  15  iterations. 
This  was  not  a  converged  design,  as  was  shown  by  the  duality  gap  between  the 
above  primal  weight  and  a  dual  weight  of  1200.66  lb.  The  optimization  may 
have  been  terminated  either  due  to  a  limit  on  number  of  iterations  or  because 
there  is  practically  no  weight  change  between  the  last  two  iterations.  The 
final  design  presented  is  close  to  the  optimum  but  has  a  significant  number 
of  constraints  (particularly  minimum  areas)  which  are  not  considered  active. 

That  this  situation  should  occur  is  puzzling  since  a  comparison 
of  the  two  methods  shows  very  few  real  differences.  In  each  case,  an  LPT 
step  (or  steps)  is  used  to  determine  a  constraint  population  followed  by  a 
N-R  search.  The  SCICON  program  allows  both  accumulation  and  discard 
of  constraints.  This  strategy  was  not  found  necessary  for  the  present  problems 
but  is  used  in  later  versions  of  the  program.  The  statement  is  made  by  RAE, 
that  the  number  of  active  constraints  is  typically  much  smaller  than  the 
number  of  variables  -  a  fact  not  really  borne  out  in  the  current  examples. 

The  RAE  solution  requires  the  solution  of  an  m  x  m  matrix,  where  m  is  the 
number  of  active  constraints.  The  m  x  m  matrix  is  itself  a  triple  product 

- 1  T 

(G  H  G  )  where  the  H  is  a  diagonal  n  x  n  matrix  with  n  as  the  number  of 
variables.  Finally,  it  should  be  noted  that  G  is  a  constraint  derivative 
matrix.  As  is  well  known,  if  the  constraints  are  stresses,  treating  the 
internal  forces  as  invariants  can  lead  to  gross  errors.  To  include  the 
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internal  force  redistribution  effects  in  computing  the  derivatives  either 
requires  a  finite  difference  method  or  the  =  -S  K  A)  formulation. 

Both  of  these  methods  are  computationally  expensive.  Due  to  the  known 
extreme  nonlinearity  of  the  problem,  it  is  entirely  possible  that  the 
convergence  characteristics  are  strongly  influenced  by  the  accuracy  of  the 
constraint  derivative  matrix  G. 


With  the  force  method,  the  constraint  derivative  matrix  is 
always  exact  and  is  generated  with  a  minimum  of  auxiliary  calculation. 


2. 8  Displacement  Constraints 

The  next  stage  is  the  formal  incorporation  of  displacement  constraints. 

In  principle,  this  does  not  raise  any  new  problems,  since  displacement 
constraints  have  already  been  treated  in  the  form  of  the  compatibility  conditions. 

The  only  obvious  feature  of  external  displacement  constraints  are 
that  they  (a)  are  associated  with  non-zero  constraint  values  and  (b)  are 
inequality  constraints  which  need  not  be  satisfied  in  the  equality  sense. 
Compatibility  conditions  are  zero-valued  and  must  be  satisfied  as  equalities 
throughout. 

The  governing  Lagrangian  is  modified  as  shown  in  Figure  9. 

Figure  9  also  gives  the  modified  forms  of  the  derivatives  (optimality 

criteria  and  constraints).  In  the  modified  equations,  one  new  matrix 
parameter  is  introduced  (Bq).  This  is  the  virtual  force  in  each  member 
of  the  structure  arising  from  the  (virtual)  unit  loads  imposed  on  the  structure 
to  calculate  the  displacements.  Strictly  from  the  definition  of  the  virtual 
load  method  Bq  need  only  be  a  system  in  static  equilibrium  with  the  virtual 
load.  In  the  present  context  it  has  been  shown'  ,  that  the  Bq  system  must 
be  the  actual  internal  force  distribution  arising  from  the  virtual  force. 

Because  of  the  additional  displacement  related  terms,  the  linearized 
Newton- Raphson  equations  must  also  be  expanded  as  indicated  in  Figure  10. 
Naturally  only  active  displacement  constraints  will  be  retained  in  the  N-R 
equations  to  be  solved.  The  total  number  of  equations  is  given  by 
(2n  +  m  +  s  -  a  +  d)  where  d  is  the  number  of  active  displacement  constraints. 
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FIGURE  10  -  MODIFIED  NEWTON  EQUATIONS  TO 

INCLUDE  DISPLACEMENT  CONSTRAINTS 


The  additional  term,  d,  docs  not  necessarily  mean  an  increase  in  the  order 
of  the  equations  since  the  total  number  of  active  constraints  can  still  never 
exceed  the  number  of  variables.  The  upper  limit  on  the  number  of  equations 
is  (2n  +  2m  -  1)  which  would  only  occur  when  no  minimum  size  constraints 
existed. 

Because  of  the  similarity  between  the  compatibility  and  displacement 
constraints,  modification  of  the  computer  program  was  relatively  straight¬ 
forward.  The  Bq  matrices  necessitated  the  introduction  of  e  ditional 
loading  cases,  but  this  concept  had  been  used  in  the  earlier  OPT1M  programs 

The  new  program  was  tested  on  small  scale  3-bar  truss  problems  as 
usual.  Figure  11  provides  details.  A  few  minor  developmental 
problems  were  encountered.  More  limits  were  introduced  on  the  early 
iterations  to  prevent  initial  transients  producing  unstable  situations.  Also 
the  ability  to  acquire  and  discard  constraints  was  also  incorporated. 

After  each  iteration,  following  the  initial  transients --estimated  as 
dying  out  in  2  or  3  iterations,  the  generated  design  was  compared  with  all 
constraints.  If  any  additional  constraints  had  become  activated,  they  were 
duly  recognized  and  considered  in  subsequent  steps.  If  the  number  of 
recognized  constraints  equals  or  exceeds  the  number  of  variables,  the 
program  would  then  return  to  theLPTmode  to  select  a  new  design.  Discard 
of  constraints  was  based  on  the  same  criterion  as  used  initially- -i.  e. ,  non¬ 
positive  Lagrangian  multipliers. 
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The  next  problem  tested  is  another  variant  of  the  10-bar  truss 


structure.  This  problem  used  in  Reference  1,  proved  to  be  a 
rather  difficult  one  to  solve  and  no  exact  optimum  solution  has  been 
published.  The  problem  is  shown  in  Figure  12.  The  major 
difficulty  arises  from  the  specification  of  constraints  on  the  vertical 
displacements  at  both  the  lower  and  upper  nodes  at  the  free  end.  If  it  is 
assumed,  say,  that  the  lower  node  will  deflect  more  than  the  upper  and 
the  structure  is  optimized  using  FSD  and  the  envelope  method  with  the 
single  displacement  constraint,  the  upper  node  in  the  optimized  system 
experiences  an  excessive  displacement.  Using  the  upper  displacement  as 
the  constrained  value  reverses  the  situation  and  the  lower  displacement 
violates  at  the  new  optimum.  Using  both  upper  and  lower  displacements 
in  the  envelope  method  produces  an  unsatisfactory  situation  in  which  no 
converged  solution  is  obtained.  In  the  course  of  the  iterations,  which  tend 
to  become  unstable  at  one  stage,  a  very  low  weight  design  is  generated. 
From  examination  of  this  design,  it  does  appear  to  be  near  optimum. 
Unfortunately  it  has  not  been  possible  to  determine  exactly  what  the  true 
optimum  is  and  which  constraints  are  active. 

Using  the  current  version  of  the  program  the  above  problem  was 
tried.  A  number  of  difficulties  were  encountered. 

In  the  first  place,  it  proved  to  be  impossible  to  generate  a  feasible 
full  vertex  design,  based  upon  any  population  predicted  by  an  LPT  solution. 
Although  a  number  of  iterations  were  permitted,  resulting  in  a  number  of 
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different  LPT  populations,  it  was  never  possible  to  converge  on  any  full- 
vertex  designs  using  the  Newton-Raphson  technique.  This  situation  had 
been  encountered  previously,  even  for  stress  limited  problems.  The 
solution  adopted  at  that  time  was  to  generate  a  fully  stressed  design  (or  a 
close  approximation  thereto)  using  a  limited  number  of  stress-ratio 
iterations.  For  this  design,  the  Lagrangian  multipliers  were  found  and  after 
any  constraint  discard,  the  full  Newton-Raphson  search  was  used.  In  the 
presence  of  potentially  active  displacement  constraints,  the  FSD  alone  may 
result  in  a  violated  design.  From  expediency,  it  was  decided  to  use  this 
FSD  suitably  scaled  to  ensure  satisfaction  of  the  single  dominant  displacement 
constraint  as  the  starting  point  for  the  Newton  search.  Consideration  was 
also  given  to  the  use  of  an  envelope  method  to  generate  a  starting  point, 
either  using  a  single  or  multiple  displacement  constraints.  This  approach 
was  not  implemented  in  the  present  program  for  a  number  of  reasons. 

Firstly,  the  envelope  method  is  computationally  more  complex  than  the 
stress-ratio  FSD.  Secondly,  from  the  previous  history  with  this  particular 
problem,  use  of  either  of  the  envelope  methods  could  be  unreliable --a  single 
constraint  always  failed  and  double  constraints  would  not  converge. 

Using  this  scaled  FSD  as  a  starting  point,  the  Newton-Raphson  stage 
was  entered.  For  a  limited  number  of  steps,  convergence  appeared  to  be 
satisfactory--moving  in  the  direction  of  the  assumed  optimal  design.  But 
after  additional  constraints  were  captured,  according  to  the  built-in  logic, 
the  direction  of  convergence  altered  and  rapidly  became  unstable. 

Table  6  provides  details. 


10-BAR  TRUSS  WITH  DISPLACEMENT  CONSTRAINTS  -  COMPUTER  RESULTS 
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After  considerable  numerical  experimentation  the  problem  appeared 
to  reside  in  the  accumulation  of  undesirable  constraints  during  the  search, 


with  no  satisfactory  indications  which  should  be  rejected.  Thus  the  problem 
again  appeared  to  devolve  on  the  question  of  constraint  identification. 
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3.0  PROGRAM  DEVELOPMENT  AND  ORGANIZATION 

The  Force  Method  Optimization  (FMO)  program  was  developed  and 
tested  using  the  Bell  RAX  system.  The  convenience  and  ready  access  of 
the  RAX  system  make  it  attractive  for  use  in  structural  and  other  program 
development,  although  the  RAX  core  storage  capacity  has  limited  the  FMO 
test  problem  size  to  no  more  than  that  required  for  solution  of  the  22-bar 
truss.  Very  recently,  a  Virtual  Storage  system  (VSPC)  has  been  installed 
at  Bell  which  considerably  increases  the  RAX  storage  capacity,  but  was  not 
yet  in  place  during  the  program  development. 

A  structure  cutter  was  also  developed  and  tested  as  part  of  the  FMO 
program  development,  but  was  not  integrated  in  the  program  since  it  would 
have  further  reduced  the  problem  size  which  could  be  analyzed  on  the  RAX 
system.  The  essential  or  basic  function  of  the  structure  cutter  in  this 
application,  i.e.,  provide  values  of  the  bQ  and  b^  matrices,  was  performed 
by  hand  calculation  and  input  as  basic  data  for  most  of  the  test  problems 
analyzed  during  program  development. 

Much  discussion  of  how  the  FMO  program  evolved  to  its  present 
form  is  sprinkled  throughout  the  Technical  Discussion,  Section  2.0,  and  the 
reader  is  referred  to  this  discussion  for  clarification  of  specific  points 
regarding  why  a  particular  action  or  step  was  required.  Also,  since  the 
FMO  RAX  program  is  reasonably  straightforward  in  programming  style, 
it  does  not  appear  to  be  necessary  to  describe  in  detail  the  programming 
steps  which  were  followed  to  implement  the  stages  of  program  development. 


Hence,  the  principal  item  to  be  discussed  in  this  section  is  the  organization 
of  the  FMO  RAX  program.  A  flow  diagram  of  the  RAX  program  is  presented 


in  Figure  13.  Referring  to  the  flow  diagram,  material  properties  and  other 
data  pertaining  to  the  structure  to  be  analyzed  are  input  to  the  program  as 
basic  data.  The  freefield  format  is  used,  as  shown  in  Figure  14,  for  the 
3-bar  truss  of  Figure  11.  The  data  to  be  entered  consists  of  the  number 
of  members,  redundancies,  and  displacement  constraints  in  the  first  line, 
followed  by  the  densities,  elastic  moduli,  and  member  lengths  in  the  second, 
third,  and  fourth  lines.  The  selected  valves  of  minimum  member  area  and 
the  upper  and  lower  stress  limits  are  entered  in  the  fifth,  sixth,  and  seventh 
lines,  followed  by  inputs  for  the  displacement  limits,  and  the  calculated 
values  of  bQ  and  bj.  This  completes  the  basic  data  entries. 

A  starting  point  design  (guess  vectors)  can  be  input  after  the  basic 
data  entry,  or  the  default  option  exercised,  which  automatically  selects 
minimum  areas  as  the  initial  point  design.  Based  on  experience,  a  closer 
approximation  to  the  optimal  system  may  be  obtained  using  a  stress-ratio 
design.  Hence,  a  stress-ratio  cycle  is  available  to  generate  a  "better" 
initial  point  design,  as  an  option. 

The  starting  or  initial  point  design  is  applied  to  the  LPTroutine  which 
generates  a  potential  constraint  population  by  solving  a  set  of  non-square 
linear  equations  in  (m  +  p)  unknowns  corresponding  to  the  optimality  criteria. 
The  LPTroutine  was  not  developed  during  this  program.  It  is  a  mathematical 
routine  obtained  from  the  Bell  library  of  computational  routines.  It  also 
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FIGURE  13.  FMO  COMPUTER  PROGRAM  -  FLOW  DIAGRAM 
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Line 


1 

3,  1,2 

No.  of  Members,  Redundancies, 
Displ.  Const. 

2 

.  1,  .  1, .  1 

Density,  lbs/cu.  ft. 

3 

10000. , 10000. ,  10000. 

Elastic  Moduli,  psi 

4 

125. , 100., 125. 

Length,  inches 

5 

.01, .01, .01 

Minimum  Area,  sq.  in. 

7 

50.  ,  50.  ,  50. 

Upper  Stress  Limit,  psi 

8 

-50. , -50. , -50. 

Lower  Stress  Limit,  psi 

9 

1. ,  .  1,  .  1 

Displacement  Limit,  in. 

10 

20. ,  -7. ,  0. 

bQ  matrix 

11 

1.  66667,  -1. 3333,  0. 

»> 

matrix 

12 

13 

0.  ,  1. , 0. 

1.  ,  -1.6,  1. 

5 


\ 

i 

s 

[i 

i 


l 


FIGURE  14  -  BASIC  DATA  -  INPUT  FORMAT 
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generates  values  for  the  non-zero  Lagrangian  multipliers.  If  the  population 
shows  that  each  member  is  constrained,  i.e.  ,  either  a  minimum  area  or  a 
maximum  stress  is  present  in  each  member,  the  solution  of  linear  equations 
leads  to  the  required  internal  force  distribution  which  satisfies  the  prescribed 
FSD  conditions  without  iteration.  If  the  constraint  population  is  other  than 
a  FSD  type,  i.e.,  a  member  may  be  defined  as  both  fully  stressed  and  at 
minimum  area,  this  will  result  in  another  member  being  completely 
unrestrained.  The  set  of  compatibility  equations  must  be  expanded  by 
additional  relationships.  The  number  of  variables  increases  to  include 
redundancies  and  undefined  areas  and  the  equations  become  nonlinear.  A 
Newton- Raphson  routine  is  used  to  generate  a  design  corresponding  to  the 
specified  constraints.  The  constraint  population  which  results  from  this 
design  is  checked  for  correspondence  with  the  original  LPT  population  output 
(optimality  criteria).  If  this  population  check  fails,  the  constrained  design 
is  iterated  (no  more  than  3  times)  through  the  LPT  routine  until  population 
agreement  is  obtained.  A  FSD  is  generated  if  population  agreement  cannot 
be  obtained  after  three  iterations,  and  the  program  passes  out  of  the  LPT  loop. 

As  noted  in  the  flow  diagram,  a  check  of  the  Lagrangian  multipliers 
is  performed  at  this  point.  If,  as  discussed  in  Sections  2.  5  and  2.  6,  all 
multipliers  associated  with  non-zero  type  constraints  have  positive  values, 
the  design  is  considered  to  be  optimal,  and  the  results  are  printed,  ending 
the  program  operation.  The  alternate  condition,  i.e.,  one  or  more  negative 
multipliers  indicates  that  non-optimal  constraints  are  present.  Discarding 
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the  constraints  associated  with  these  negative  multipliers  leaves  a  reduced 
set  of  potential  active  constraints  for  the  optimum  design.  The  resulting 
nonlinear  equations  are  solved  iteratively  using  the  Newton- Raphs on 
routine  until  convergence  is  obtained,  indicating  an  optimum  design  has 
been  reached.  Results  are  printed,  ending  program  operation. 


CONCLUSIONS 


4.  0 

The  new  approach  to  structural  optimization  via  the  use  of  the  force 
method  of  analysis  and  a  linear  programming  stage  has  achieved  a  major 
degree  of  success.  Some  of  the  classic  optimization  test  problems  have 
been  reduced  to  almost  trivial  solutions.  The  10-bar  truss  problem  with 
non-uniform  stresses  is  solved  without  iteration.  On  the  larger  22-bar 
trusses,  solutions  have  been  achieved  in  the  majority  of  cases,  which  are 
confirmed  as  being  optimal  by  comparison  with  analytically  generated 
results.  The  results  provide  a  clear  indication  of  the  validity  of  the  approach 
for  the  rigorous  incorporation  of  stress  constraints.  Because  of  the  relatively 
rapid  convergences  achieved,  the  iterative  computational  effort  appears  to 
be  significantly  lower  than  that  required  by  alternate  approaches.  The  key 
to  the  generation  of  the  optimal  design  lies  in  the  ability  to  predict  or 
identify  the  active  constraints.  The  linear  programming  stage  incorporated 
in  the  developed  approach  does  an  excellent  job  of  selecting  dominant 
constraints  through  local  linearization  of  the  non-linear  domain,  but,  as  in 
any  non-linear  situation,  does  require  either  a  good  starting  point  and/or 
a  well  conditioned  problem  to  solve. 

In  the  study  of  the  22-bar  problems,  no  solutions  could  be  generated 
using  the  standard  program,  although  the  correctness  of  the  optima  were 
demonstrated  by  the  program.  In  these  cases,  one  of  which  was  a  full 
vertex  design,  the  problem  was  felt  to  be  attributable  to  conditioning.  This 
could  not  be  entirely  proven,  since  some  of  the  other  22-bar  cases  which 
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were  successfully  solved  might  have  appeared  to  be  more  poorly  conditioned. 
The  conditioning  problem  here  was  taken  to  be  associated  with  the  very  large 


spreads  in  magnitude  between  the  largest  and  smallest  Lagrangian  multipliers. 
The  spread  is  10^  and  since  the  merit  function  is  the  sum  of  the  multipliers, 
there  is  potential  for  round-off  error  to  obscure  the  situation  and  lead  to 
erroneous  definition  of  active  constraints. 

Naturally,  the  22 -bar  problem  is  an  unrealistic  one  selected  to 
exhibit  certain  characteristics  on  a  relatively  small  scale  and  to  be  a  test 
for  optimization  techniques.  It  is  not  clear  whether  the  problem  is  really 
pathological  -  especially  in  its  ranges  of  values  of  multipliers  -  or  if  it  is 
truly  characteristic  of  larger  scale  problems.  Some  consideration  has  been 
given  to  transformation  of  the  problem  to  try  and  eliminate  the  size  disparity 
in  numbers  which  are  linearly  combined.  Simple  scaling  alone  will  not 
suffice,  but  no  suitable  transformation  could  be  evolved. 

The  presence  of  very  small,  but  non-zero,  values  for  Lagrangian 
multipliers  associated  with  real  (active)  constraints  presents  a  hazard  in 
the  Newton-Raphson  search  procedure  which  allows  for  constraint  acquisition 
and  discard.  The  question  of  when  a  multiplier  becomes  zero  and  the  corres¬ 
ponding  constraint  should  be  discarded  is  made  extremely  difficult,  since 
such  small  values  are  readily  masked  by  round-off  error  in  large  scale 
problems. 

With  regard  to  displacement  constraints,  the  fundamental  concepts 
and  approach  seem  to  have  been  validated.  Some  small  scale  problems  were 
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successfully  solved.  Difficulties  with  constraint  detection  on  larger 
problems  appear  to  remain.  In  purely  stress  and  minimum  size  problems, 
the  number  of  active  constraints  equals  or  approaches  the  number  of 
variables.  The  population  defined  by  the  LPT  stage  is  therefore  a  reliable 
indication  of  the  final  active  constraints.  In  the  presence  of  displacement 
constraints,  this  no  longer  holds  true.  The  number  of  active  constraints 
is  frequently  very  much  reduced  and  the  LPT  is  much  less  reliable  in 
predicting  population. 

To  overcome  this  situation,  the  expedient  solution  of  selecting  a 
single  dominant  displacement  constraint  was  adopted  to  generate  a  starting 
design  for  the  Newton  search.  This  did  not  perform  satisfactorily  and 
placed  too  much  reliance  on  the  constraint  acquisition/discard  algorithms. 
It  is  clear  that  the  method  of  selecting  a  starting  design  for  the  Newton 
search  must  be  revised,  but  in  what  direction  requires  some  further 
consideration. 

Although  difficulties  were  encountered  in  the  incorporation  of 
displacement  constraints  per  se,  it  is  believed  that  these  are  operational 
problems  of  technique  and  strategy  rather  than  of  concept. 

In  all  optimization  research,  there  has  always  been  a  great  deal 
of  concern  for  computational  efforts  involved  in  the  various  approaches. 
Both  the  numbers  of  iterations  and  the  cost  per  iteration  are  major 
contributing  factors.  In  this  force  method  work,  the  costs  appear  to 
compare  most  favorably  with  other  methods.  The  force  method,  generally, 


requires  the  solution  of  a  smaller  array  of  equations  than  the  displacement 
method  for  a  given  structure  (redundancies  versus  displacements).  The  need 
for  a  structure  cutter  in  the  force  method  has  a  major  cost  impact  in  a 
single -shot  analysis  but  this  effect  is  reduced  considerably  in  iterative 
situations.  The  LPT  stage  is  equivalent  to  the  inversion  of  an  array  of 
order  of  the  number  of  variables.  The  full  Newton-Raphson  search  does 
require  consideration  of  a  larger  set  of  equations  -  iteratively.  For  the 
22-bar  trusses  31  equations  are  required.  Although  this  array  is  large, 
the  computation  of  the  individual  terms  is  very  straightforward.  No 
auxiliary  analyses  are  required  to  generate  approximations  to  gradients 
of  constraints.  The  absence  of  approximations  throughout  is  considered 
to  be  a  major  factor  in  the  good  convergence  characteristics  observed. 

The  size  of  the  Newton  array  may  be  a  problem  in  larger  systems.  The 
matrix  is  symmetric  and  sparse,  so  that  the  possibility  of  using  efficient 
solution  techniques  exists. 

In  view  of  the  successful  solutions  demonstrated,  it  is  recommended 
that  development  of  this  force  method  approach  be  continued.  Clearly,  the 
problems  encountered  will  have  to  be  overcome.  Further,  study  of  the 
conditioning  problem  in  the  LPT  stage  may  be  solvable  either  through  some 
transformation  or  use  of  a  more  sophisticated  LPT  routine. 

For  displacement  constraints,  a  method  of  selecting  a  better  starting 
point  for  the  N.  R.  search  is  required  when  the  LPT  is  unsatisfactory.  Also 
the  strategies  for  acquisition  and  discard  of  constraints  during  the  N.  R. 
stage  may  require  further  treatment. 
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Finally,  the  search  for  alternate  methods  of  handling  the  entire 
problem  of  constraint  identification  should  be  pursued.  The  present  work 
is  felt  to  have  been  a  positive  step  in  the  direction  of  generating  an 
operational  optimization  program  and  provides  a  firm  basis  for  future 
work  in  this  complex  but  rewarding  technology. 
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